{
 "metadata": {
  "name": "testing"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import powerlaw\n",
      "%load_ext rmagic "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "During startup - Warning messages:\n",
        "1: package \u2018methods\u2019 was built under R version 2.15.1 \n",
        "2: package \u2018datasets\u2019 was built under R version 2.15.1 \n",
        "3: package \u2018utils\u2019 was built under R version 2.15.1 \n",
        "4: package \u2018grDevices\u2019 was built under R version 2.15.1 \n",
        "5: package \u2018graphics\u2019 was built under R version 2.15.1 \n",
        "6: package \u2018stats\u2019 was built under R version 2.15.1 \n"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "data = genfromtxt('words.txt')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 14
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "data = genfromtxt('words.txt')\n",
      "####\n",
      "fit = powerlaw.Fit(data, discrete=True, estimate_discrete=False)\n",
      "print fit.xmin\n",
      "print fit.alpha\n",
      "print fit.power_law._pdf_discrete_normalizer\n",
      "print fit.distribution_compare('power_law', 'lognormal')\n",
      "\n",
      "data = genfromtxt('words.txt')*10.0**2\n",
      "####\n",
      "fit = powerlaw.Fit(data, discrete=True, xmin=700.0, estimate_discrete=False)\n",
      "print fit.xmin\n",
      "print fit.alpha\n",
      "print fit.power_law._pdf_discrete_normalizer\n",
      "print fit.distribution_compare('power_law', 'lognormal')\n",
      "\n",
      "data = genfromtxt('words.txt')*10.0**1\n",
      "####\n",
      "fit = powerlaw.Fit(data, discrete=True, xmin=70.0, estimate_discrete=False)\n",
      "print fit.xmin\n",
      "print fit.alpha\n",
      "print fit.power_law._pdf_discrete_normalizer\n",
      "print fit.distribution_compare('power_law', 'lognormal')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Calculating best minimal value for power law fit\n",
        "7.0"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "1.95271772619\n",
        "5.67848426332"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "(0.93878265949091055, 0.42433297291177607)"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "700.0"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "2.02133980333\n",
        "821.609142112"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "(9.5167210183711806, 2.3670503345783636e-07)"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "70.0"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "2.01472443762\n",
        "75.069370313"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "(6.5019676042704599, 0.00046635463846706013)"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      }
     ],
     "prompt_number": 48
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.special import zeta\n",
      "zeta(2.0, 1.0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 56,
       "text": [
        "1.6449340668482266"
       ]
      }
     ],
     "prompt_number": 56
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "zeta(2.0, 10.0)/zeta(2.0, 1.0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 69,
       "text": [
        "0.063933465663574926"
       ]
      }
     ],
     "prompt_number": 69
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "zeta(2.0, 1.0) - zeta(2.0, 10.0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 63,
       "text": [
        "1.539767731166541"
       ]
      }
     ],
     "prompt_number": 63
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = 0\n",
      "for i in range(1, 10):\n",
      "    x += i**-2 \n",
      "print x"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1.53976773117\n"
       ]
      }
     ],
     "prompt_number": 66
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from numpy import log, sum\n",
      "data = genfromtxt('words.txt')\n",
      "factor = 100.0\n",
      "xmin = factor*7.0\n",
      "data = data*factor\n",
      "data = data[data>=xmin]\n",
      "n = len(data)\n",
      "print 1 + ( n / sum( log( data / ( xmin - 0.5) ) ))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "2.02138372655\n"
       ]
      }
     ],
     "prompt_number": 46
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%prun\n",
      "data = genfromtxt('cities.txt')/10**3\n",
      "####\n",
      "import powerlaw\n",
      "results = powerlaw.Fit(data, discrete=False)\n",
      "results.xmin\n",
      "results.fixed_xmin\n",
      "results.alpha\n",
      "results.D\n",
      "len(results.xmins)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Calculating best minimal value for power law fit\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#data = genfromtxt('flares.txt')\n",
      "#discrete=False\n",
      "data = genfromtxt('quakes.txt')\n",
      "data = (10.0**data)/(10.0**3)\n",
      "discrete=False\n",
      "#data = genfromtxt('terrorism.txt')\n",
      "#discrete=True\n",
      "#data = genfromtxt('words.txt')\n",
      "#discrete=True"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 33
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fit = powerlaw.Fit(data, discrete=discrete, discrete_approximation='round')#15000)\n",
      "#data = fit.data\n",
      "xmin = fit.xmin\n",
      "print fit.xmin\n",
      "print fit.alpha"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Calculating best minimal value for power law fit\n",
        "10.0"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "1.95240332155\n"
       ]
      }
     ],
     "prompt_number": 34
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fit.plot_ccdf(original_data=True)\n",
      "fit.D"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 44,
       "text": [
        "0.038097679127198414"
       ]
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEHCAYAAABGNUbLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHzpJREFUeJzt3XtclVW+x/EP6HiLohyV9GhUhokTFlNeC9iNBPQiprIs\nOWewwBS1I6ZlntFKyl41Nk6pjKXlZbrpdCo7WqYFOWzQknTUIsvo4q5Ja8pMahxDjef8sfKCom42\nz97Pvnzfr9fzQh7Z+/mK8Ntrr7WetaIsy7IQEZGIEu10ABERCTwVfxGRCKTiLyISgVT8RUQikIq/\niEgEUvEXEYlAKv4iIhFIxV9EJAK19NcT19XV8fvf/569e/dy9dVXk5WV5a9LiYhIE/mt5b927Vr6\n9OnDY489xtKlS/11GRER8UGTin9BQQFxcXEkJSU1OF9RUUFiYiIJCQmUlJQAUF1dTffu3QHYu3ev\nTXFFRMQOTSr++fn5rFq16pjz48aNY968eZSVlTFnzhx27txJ7969+fTTTwFo166dPWlFRMQWTerz\nT0lJwePxNDhXW1sLQGpqKgAZGRlUVVWRnp7OlClTWLt2LYMHD270+aKionyILCIizV2Ts9l9/uvX\nr6dnz56HPu/Vqxfr1q2jdevWzJgxg5KSEjIzM4/7+NWrLRYvtnj4YYs777QYNswiI8Oid2+LTp0s\nWrY0H3v3NucvvHAqd95pvn7xYovVqy3ef99i1y6L+noLyzr2mDp1aqPnj/d3R5870eOb8rzNeUyw\n5DzZ10dyzkD8nyunvTlD9XfdDn6b7eMtt7sYl8tFbq6r0b8/cAB27oSvvjJHebmL00+Hzz6DqqrD\n57/6Cn78EeLi4Je/hBYtDj/HDz+4WLGi8es39ncHzx18Y7J9O7z6qvnzwXNRUYePg59HRx/+WFvr\norLy8LmD51u0MB+PPFq0MMfOnS7y8+EXvzh8tGoFrVvD9u0uHnoI2rQxn7dpA5blYtky8+c2bUzO\n996DM84wR9u2h/Mdj8vV+Pfd169v7O+beg1frtvUr/dHTl8ef6LHHO/vlNO3xwTLz+bJHnOinOXl\n5ZSXlzf5eo2Jspr4MuLxeMjJyaG6uhow3T4ul4tNmzYBMHbsWLKyssjOzj75xaOibHsVA1P8//lP\n+OYbaO7THvn4J54oZsSI4kPnLOvwcfTn9fXm408/HXuuvr7h8dNPDT8eOGCO/fsPH/v2QV3d4Y91\ndebfWVcHe/ce/vzHH+GTT4qJjS3mu+/gu+/Mc55xBpx+esOjfXvzAnnwY8eO0KmTOTp2NC8k/lRc\nXExxcbF/L2ID5bRXKOQMhYxgT+1sdss/NjYWMDN+zjrrLEpLS5k6darXjy8uNi1/O16B27SB+Hhz\n2Onf/3bRt6+9z+kP5eUujvw2/vijeRHYvfvw8d13sGuXOT77DDZuNC+WX399+GjXDrp0gc6dzccL\nLoDBgyEhwZ6cdvxfB4Jy2isUcgZ7Rsda/rm5ubjdbr799ls6derEfffdR35+Pm63m1GjRrF//36K\nioooKiry7uI2t/yl+SzLvEDs2AFffmm6kt5+G156CTp0gOuuMy8ESUkn71ISEf+wo3Y2udvHTir+\noaO+Ht56C5YuNUd0NFxzDVx7LQwY0HCMRUT8y47a6fjaPsXFxba9jRH/iY6GSy+FP/0JPv0UXnwR\nTj0Vbr3VdA2NHAmvvWbGJkTEP8rLy20bk1DLX5rtk09Mt9CLL8KHH8JVV0FuLlxxBbR0fD6ZSPhR\nt48EnS++MN1Czz5rBpRzcyEvD5KTNUYgYhd1+0jQ6doViorMPRhut+kauv56M2PowQfNC4KI+Ebd\nPhJSLAvefBOefhqef968EEyYAFdf7XQykdCkbh8JOXV1sGIFTJ4Mv/oVlJSYAWMR8Z66fSTktG5t\n7hPYvBkSE+HCC+GJJ5p/R7ZIJFC3j4SNd9+FW24x6xg9+CD8vDisiJxAWLT8JbL17m1uHisshJtv\nhsxMWL/e6VQi4U/FXxzXogUMGwZbt5o7hg8eH3zgdDKR8KXiL0GjVSsYNQo++sjcTZyaat4RfPml\n08lEwo/jxV8DvnK0tm3hjjvM3cKnnWamht59N3z/vdPJRJylAV+JKJ99Zor/66+bjyNHmgFikUil\nAV+JCPHx8NRTsGoVLF9u7g944QVNDxVpDrX8JeS8/jrceafZdGbGDBg40OlEIoGllr9EpIwMswPZ\nqFEwdKhZO+ijj5xOJRJaHC/+GvAVX0RHm+mhH34IF19sNpS5/351BUl404CvyFF27DA7i/XoAfPn\n+38TehEnqdtH5GddukB5uVk4btAgsxG9iByfir+EjXbt4Lnn4De/gf794b33nE4kErxU/CWsREfD\ntGlw773gcsGf/6xxAJHGqM9fwlZNDfzud9C+PSxcqH0DJHyoz1/kBHr0gLVrzUyg5GSzi5iIGGr5\nS0R4+23zLmDAAJg9G2JjnU4k4ruwaPlrnr8EQt++sGmTWTTuoougstLpRCJNp3n+Is3w8stmcbib\nb4b77tMicRJ6tIG7iI++/hpuugn+9S8zPVSDwRJKwqLbR8QJnTrBihVmnaBLLjE3iIlEErX8JeK9\n/rpZJ2j8eJg40dwrIBLM1O0jYpN//ANuvNHcJbxoEXTr5nQikeNTt4+ITbp1g4oKuPxys0roU0/p\nzmAJb2r5ixxl0ybTDXTeebBggblDWCSYqOUv4gfJybBhA5x5Jlx9tVkpVCTc+LX4b9u2jVtuuYUh\nQ4b48zIitmvdGubMgY4dYfRodQFJ+PFr8T/nnHOYP3++Py8h4jfR0abvf9Mm+NOfnE4jYi+vin9B\nQQFxcXEkJSU1OF9RUUFiYiIJCQmUlJT4JaCIk2JiYPlyeOQReOUVp9OI2Mer4p+fn8+qVauOOT9u\n3DjmzZtHWVkZc+bMYefOnTz99NOMHz+eHTt22B5WxAndusHSpZCfD+++63QaEXu09OaLUlJS8Hg8\nDc7V1tYCkJqaCkBGRgZVVVXk5eWRl5cHwK5du5g8eTKbN29m+vTpTJo06ZjnPnKRIpfLhcvl8uGf\nIeJf/frBo4/CFVfAM8+YjyKBUl5ebvsCmF5P9fR4POTk5FBdXQ1AWVkZCxYsYMmSJQDMnTuX7du3\nM23aNO8vrqmeEmLcbnMz2N13w623Op1GIlVYTPXUks4SStLS4M03zbuAW2+F/fudTiSRxJElnY9u\n+dfW1uJyudi0aRMAY8eOJSsri+zsbO8vrpa/hKjvv4ehQ03xf/55OP10pxNJJHG05R/781ZIFRUV\neDweSktL6devX7PCiISK004z+wL86lcwcCBs2+Z0IpGm8ar45+bmMnDgQGpqaujWrRuLFi0CYObM\nmRQWFpKens6YMWPo0KFDkwOo20dCVYsWMHMmjBljXgDeesvpRBLutJOXSJBZscJMBS0pMQPCIv5k\nR+30aqqnPxUXF2uKp4S87GwoLYWcHPj4Y5g8GaKinE4l4cbOKZ9q+YvYaMcO8wKQlASPPw6tWjmd\nSMJRWEz1FAknXbqYfQF27zZbRO7a5XQikcY5Xvw14Cvh5pRT4MUXzaYwAwbARx85nUjChQZ8RULE\nvHlwzz3wv/9rbhATsYO6fUSCXGEhPPssDBli9gYWCRZq+YsEwAcfwFVXmReBBx4wewWI+CosWv7q\n85dIkJgIVVVQWQlFRdoZTHyjPn+REFVba/r+r7vOrAwq4ouwuMlLJJLExsKqVXDppdCpkxkTEHGC\nir9IgJ15Jrz2GqSmmg3iBw92OpFEIseLv5Z3kEh03nlmT+DMTHNfQGam04kkFGh5B5EwUVlpZgBN\nnAgTJmg9IPGOHbVTxV/EYZ99Zrp+zj8f5s+Hdu2cTiTBLiymeopEuvh4WLMGfvELsxzEp586nUgi\ngYq/SBBo2xb+8hcYPhwuuww2b3Y6kYQ7dfuIBJkXXjCbw7/0ktkhTORoYdHtozt8RRq6/np46im4\n5hp4/XWn00gw0R2+IhFg7VozEDxnjnlBEDlIs31EwtzmzWZBuMJCmDJFC8KJoeIvEgG+/NK0/Dt1\nMt1Bp57qdCJxWlj0+YvIiXXuDH/7G8TFQb9+UFPjdCIJByr+IiGgVSuYOxfGjzeLwj3yCOzf73Qq\nCWUq/iIhZMQIsyTE669D795mhVARX6jPXyQEWRa8+qp5J9CjByxcaMYEJDKERZ+/5vmLNF1UFGRn\nw3vvQc+ekJUF33/vdCrxN83zF5FDLAv++79hyxbTDdSmjdOJxN801VNEAKivh//6L9i71ywP0dLx\nnTrEn8Ki20dEmi86Gp58EurqzKBwfb3TiSTYqfiLhIlWrUyrv6bGvAvQGICciIq/SBg55RQoLTUb\nxScnw/r1TieSYKXiLxJm2rUzN4Q99JCZETRjhrqB5Fga8BUJY599Bv/5n+YFYf58s2uYhD4N+IrI\nCcXHg9sNgwbBJZfA44+bqaEifm35L1u2jBUrVnDgwAFGjRpF3759G15cLX+RgNmyBfLzzXiA3gWE\ntpCZ5//1118zdepUHnvssYYXV/EXCagDB8wYwCOPwPPPQ2qq04nEFwHr9ikoKCAuLo6kpKQG5ysq\nKkhMTCQhIYGSkpLjPn769OkUFhY2K6iINF/LlvA//wPPPmv2CPjrX51OJE7xquVfWVlJTEwMw4YN\no7q6+tD55ORkZs2aRXx8PJmZmaxZs4aVK1eyceNGJk6cSOfOnZk0aRKZmZkMGjTo2Iur5S/imOpq\nMxvo1lvhzjvNekESGuyonV7dBJ6SkoLH42lwrra2FoDUn983ZmRkUFVVRV5eHnl5eQDMnj2b1atX\n88MPP/Dxxx832vo/cpEil8uFy+Xy4Z8hIk2VlARvvWVeALZtg9mzzY1iEnzKy8ttXwDT6z5/j8dD\nTk7OoZZ/WVkZCxYsYMmSJQDMnTuX7du3M23aNO8vrpa/iOO+/x7y8uCrr+C55+Dss51OJCcTFlM9\ntaSziLNOOw3+7//ghhvMNpHLlzudSI7HkSWdj27519bW4nK52LRpEwBjx44lKyuL7Oxs7y+ulr9I\nUHnrLRg6FIYMgenToUULpxNJYxxt+cfGxgJmxo/H46G0tJR+/fo1+XnU8hcJHgMGwMaNsGEDjBun\nG8KCTcBb/rm5ubjdbr799ls6derEfffdR35+Pm63m1GjRrF//36KioooKipq2sXV8hcJSrW1cNll\n5qawCROcTiNHC5mbvI57cRV/kaD1j3+YdwKzZsF11zmdRo4UsKme/lRcXKwpniJBqFs3ePllyMyE\nLl3MC4E4y84pn2r5i8gJrVwJBQVmf+ALL3Q6jUCYTPUUkeB25ZUwcyakp8PDD2tvgHDhePHXbB+R\n4HfjjVBVBS++CFdcAV984XSiyOTIPH9/ULePSGg5cMDM/581yywHMXSo04kik2b7iIgj1q+HYcPM\n+kCPPgodOjidKLKERZ+/un1EQk+fPuZmsLPOgt69tSREoKjbR0SCRmUl3HwzpKWZrqCYGKcThb+w\naPmLSGhLSYF33jFLQVxyCbz7rtOJxBsq/iLSbDExsGgR3HWX2Sz+sce0LlCwU7ePiNiqpsZMDY2P\nhwcegF69nE4UfsKi20cDviLhpUcPszR0375w+eVmr+CfV36XZtKAr4iEhD174PHHYcYM+PWvYfFi\nOPVUp1OFPs3zF5GQ8OOPMHIkREfDX/7idJrQFxbdPiIS/tq0MTeDvfmm2SdYnKeWv4gEzN//bhaK\nW7/eDAiLb8Ki5a8BX5HIcfHFcMcd8LvfmXWCpGk04CsiIau+3qwM6nLB3Xc7nSY0acBXRELS9u3m\nXcALL5i9gqVpwqLbR0Qiz3/8h7kjeOhQ+PJLp9NEJhV/EXHElVfCiBHmbuD9+51OE3nU7SMijqmv\nh5wcSEgwW0WKd9TtIyIhLToannkGXn4ZlixxOk1kUctfRBz3zjtmg/hVq8xAsJxYWLT8Nc9fRC68\nEJ54ArKytBz0iWiev4iEpYPLQZ97LsyfD2ec4XSi4BQWLX8RkYN69IB166BrV0hOBrfb6UThSy1/\nEQlKL78Mt94K/fvDH/5g3g2IoZa/iIStnBzYutWMB/TtC3feCbt3O50qfKj4i0jQatcOpkyB996D\n776D3r3ho4+cThUe1O0jIiHjiSfgvvtg9WpzY1iksqN2trQpi4iI340YAVFR8JvfwBtvmAFi8Y2K\nv4iElFtuMS8AgwbpBaA5/Fr8t27dyqxZs9i3bx/Z2dkMHjzYn5cTkQgxfPjhdwDl5XDeeU4nCj0B\n6fPft28fN910E0uOWrxDff4i0hxPPAH332/uBzj7bKfTBE7ApnoWFBQQFxdHUlJSg/MVFRUkJiaS\nkJBASUlJo49dvnw5l19+OTfccEOzgoqIHG3ECLMt5KBB8MUXTqcJLV61/CsrK4mJiWHYsGFUV1cf\nOp+cnMysWbOIj48nMzOTNWvWsHLlSjZu3MjEiRPp0qXLoa/97W9/y/LlyxteXC1/EbHBjBnmXYDb\nDWee6XQa/wvYbJ+UlBQ8Hk+Dc7W1tQCkpqYCkJGRQVVVFXl5eeTl5QHgdrtZunQplmUxZMiQZgUV\nETmeO+6AvXvNyqCLF5v7AeTEfB7wXb9+PT179jz0ea9evVi3bh3Z2dmHzqWlpZGWlnbC5zlyhTqX\ny4XL5fI1kohEsLvugl/+0qwM2rev2Rw+XJaHLi8vt331Y8enetq1PKmIRLaoKBgzBvLzzYqg11wD\nSUnw8MNwRDs1JB3dML733nub/Zw+L+/Qp08ftm7deujzLVu20L9//yY/j9bzFxE7tW0LY8fCxx9D\nZiakpUFZmdOp7OHIev4ej4ecnJxGB3zPOusssrKyWLNmDR06dPD+4hrwFRE/c7vNHgH33guFhU6n\nsUfApnrm5uYycOBAampq6NatG4sWLQJg5syZFBYWkp6ezpgxY5pU+A9Sy19E/CktDSorTffP+PHw\n009OJ/KddvISEWmi776D664zM4FmznQ6TfPYUTtV/EUkYnzzjRn8fftt6N7d6TS+C4vNXNTtIyKB\n0rEjjBsH99zjdBLfqNtHRMRH//qXWQhu1Sq46CKn0/gmLFr+IiKBFBNjdgebPNnpJM5yvPir20dE\nAm3kSPjgAzMNNJSo20dEpJmeeQYefRTWrjV3B4cSdfuIiPgoN9f0/x+12HDEcLz4q9tHRJzQooVZ\nCrqgwIwB7NzpdKKTU7ePiIhNPB74wx/g+efNC8Httwf/ngDq9hERaaazz4a5c+Gdd6CuDi68EHbt\ncjqV/6nlLyJyhJEj4YwzYPp0p5McX1i0/NXnLyLBZOpUsx/A9u1OJzmW+vxFRPxo0iTYvRvmzXM6\nSeO0sJuIiB/s2gXnn2/uAejRw+k0xwqLbh8RkWDTvj1MmGD2BQ5XavmLiDTi3/+GhARYtgwuucTp\nNA2p5S8i4ift2sHdd4fvAnCOF3/N9hGRYDV8OHz9tdkQ/sABp9Noto+ISMDs3g033GCWg/jrXyE2\n1ulE6vYREfG700+HV1+Fc86BSy+FbducTmQPFX8RkZNo2RLmzIHCQhg4ED75xOlEzaduHxGRJhg9\n2mwDefvtzmVQt4+ISIClpUFlpdMpmk8tfxGRJvjiC7Px+zffOLcDmFr+IiIB1rWr2QT+ww+dTtI8\njhd/zfMXkVCTkgJr1gT+uprnLyLioHnz4M034cknnbm+un1ERBzgVMvfTir+IiJN1LOnufN3xw6n\nk/hOxV9EpImio83dvqHc+lfxFxHxQah3/aj4i4j44LLLQrv4a7aPiIgP9u0zO37t2AGnnRbYa2u2\nj4iIQ1q1Mjt8vfWW00l84/fiv2fPHvr06cOKFSv8fSkRkYAK5a4fvxf/hx56iBtvvNHflxERCbiw\nL/4FBQXExcWRlJTU4HxFRQWJiYkkJCRQUlJyzONKS0vp1asXHTt2tCetiEgQGTgQNmww/f+hxqsB\n38rKSmJiYhg2bBjV1dWHzicnJzNr1izi4+PJzMxkzZo1rFy5ko0bNzJx4kQeffRR9uzZw/vvv0/b\ntm156aWXiDpiGTwN+IpIqEtOhsceg/79A3dNO2pnS2++KCUlBY/H0+BcbW0tAKmpqQBkZGRQVVVF\nXl4eeXl5ANx///0APPnkk3Ts2LFB4T/oyEWKXC4XLperqf8GERHHuFzwwAOwcCF06OCfa5SXl9u+\nAKbXUz09Hg85OTmHWv5lZWUsWLCAJUuWADB37ly2b9/OtGnTvL+4Wv4iEuL27IG77oIlS+DhhyE3\n1//r/IfFVE8t6SwioeyUU+CRR2DZMnjwQbjqKvj8c/9cy84lnX0u/n369GHr1q2HPt+yZQv9fej0\nKi4uVlePiIS8fv3g73+HAQPg17+GP/8ZfvrJ3mu4XC7ni39sbCxgZvx4PB5KS0vp16+fLaFEREJR\nq1amC2jNGnjuObP+z/vvO52qcV4V/9zcXAYOHEhNTQ3dunVj0aJFAMycOZPCwkLS09MZM2YMHXwY\n7VC3j4iEm549we2GvDyz4fu990JdXfOfVzt5iYiEiC++gDFj4JNPYP580y3UXBrwFREJcl27msHg\nqVNh8GAoKoIffvDtudTyFxEJQbt2we23w+rVMHcuXHmlb89jR+1U8RcRCbCyMhg50nQBzZwJTV0B\nR90+IiIhKD0dqquhc2dISoJnngFvarm6fUREwsSGDTB8OHTpYrqC4uNP/piwaPmLiESySy4xLwCp\nqXDxxTB7tv03hzVGLX8RkSBRUwMjRph7AhYuhF69Gv+6sGj5q89fRMTo0QP+9jcoKGh8Oqj6/EVE\nIlhYtPxFRCTwVPxFRCKQ48Vfff4iIt5Rn7+ISARTn7+IiPhExV9EJAKp+IuIRCAVfxGRCOR48dds\nHxER72i2j4hIBNNsHxER8YmKv4hIBFLxFxGJQCr+IiIRSMVfRCQCqfiLiEQgx4u/5vmLiHhH8/xF\nRCKY5vmLiIhPVPxFRCKQir+ISARS8RcRiUAq/iIiEUjFX0QkAqn4i4hEIL8W//LyclJSUhg9ejRu\nt9ufl/KrULkJTTntpZz2CoWcoZDRLn4t/tHR0cTExNC6dWvOPfdcf17Kr0LlB0I57aWc9gqFnKGQ\n0S5eFf+CggLi4uJISkpqcL6iooLExEQSEhIoKSk55nEpKSmsXLmS2267jRkzZtgS2Jf/nBM9prG/\ns+MHIFxznuzrIzlnIP7Pfb1Ocx8frjkj6Xf9aF4V//z8fFatWnXM+XHjxjFv3jzKysqYM2cOO3fu\n5Omnn2b8+PHs2LGDqKgoANq3b8+ePXtsCRwq3+hwzani7/s1m/qYUC6qvl6nuY8Phd+hkz0mUMUf\ny0vbtm2zLrjggkOf796927rooosOfT527FjrlVdeafCYpUuXWoWFhdbw4cOtDRs2HPOcgA4dOnTo\n8OForpb4aP369fTs2fPQ57169WLdunVkZ2cfOnfttddy7bXXHvc5LC3qJiLiCE31FBGJQD4X/z59\n+rB169ZDn2/ZsoX+/fvbEkpERPzL5+IfGxsLmBk/Ho+H0tJS+vXrZ1swERHxH6+Kf25uLgMHDqSm\npoZu3bqxaNEiAGbOnElhYSHp6emMGTOGDh06+DWsiIjYw9GdvERExBlBOeC7detWRo8ezfDhw1m6\ndKnTcY5r2bJljBw5koKCAt5++22n4xzXtm3buOWWWxgyZIjTURpVV1fHhAkTGD16dKP3kwSLYP8+\nHhQqP5eh8nsOsGfPHvr06cOKFSucjnJc5U1dTqfZk0X9qK6uzho6dKjTMU7qn//8pzVq1CinY5zU\n9ddf73SERr3xxhvW4sWLLcuyrBEjRjic5uSC9ft4tFD5uQyF3/N77rnH+uMf/3jMvUzBxO12W1lZ\nWda4ceOszz///KRf79eWv6/LQgAsX76cyy+/nBtuuMGfEZudE2D69OkUFhb6O2azcwZSU7JWV1fT\nvXt3APbu3Ru0OZ3kS85A/Vweqak5A/l77mvO0tJSevXqRceOHQOasak5m7ycjj9fiSoqKqyNGzc2\nuDPYsizroosustxut+XxeKzzzz/f+uabb6ynnnrKuu2226zt27c3+NqcnBx/RmxWzvr6emvixIlW\nWVmZ3zM2J+dBgWyxNiXr6tWrrSVLlliWZVkjR44MWMam5jzIiZa/tzl37twZ8J9LX3Ie+f20rMD8\nnh+pKTmnTJli3XbbbVZGRoZ19dVXW/X19UGZ86Da2lpr+PDhJ31un+/w9UZKSgoej6fBudraWgBS\nU1MByMjIoKqqiry8PPLy8gBwu90sXboUy7IC0r/qa87Zs2ezevVqfvjhBz7++GO/t7J8zblr1y4m\nT57M5s2bmT59OpMmTfJrzqZmTU9PZ8qUKaxdu5bBgwf7PZuvOQcMGBDw72NTc65bt45PP/00oD+X\nvuSsqqoiJiYmoL/nvua8//77AXjyySfp2LHjoTXLgi3nvn37eO211zhw4ACjR48+6XP7tfg3xptl\nIdLS0khLSwt0tAa8yVlUVERRUZET8Q7xJmf79u2ZO3euE/EaOFFWu1Z9tcOJcgbD9/Gg4+WcNm0a\nY8eOdTBZQyfK6fTv+ZFO9rt00003ORWtgRN9P0+0nM7RgnK2j4iI+FfAi3+oLAuhnPYLlazKaS/l\ntJddOQNe/ENlWQjltF+oZFVOeymnvWzLafPgdANDhw61OnfubLVq1crq2rWrtXDhQsuyLKu8vNzq\n2bOn1b17d2vWrFn+jOAV5bRfqGRVTnspp738mVPLO4iIRCAN+IqIRCAVfxGRCKTiLyISgVT8RUQi\nkIq/iEgEUvEXEYlAKv4iIhFIxV9EJAKp+IuIRKD/B6h4XVnhfgEdAAAAAElFTkSuQmCC\n"
      }
     ],
     "prompt_number": 44
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "bins, CDF = fit.cdf()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 47
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fit.n_tail"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 49,
       "text": [
        "2659.0"
       ]
      }
     ],
     "prompt_number": 49
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "n = 38\n",
      "plot(fit.xmins[:n], fit.Ds[:n])\n",
      "plot(fit.xmins[:n], fit.sigmas[:n])\n",
      "plot((.794, .794), ylim())\n",
      "plot((10, 10), ylim())"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 39,
       "text": [
        "[<matplotlib.lines.Line2D at 0x82b75d0>]"
       ]
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAD9CAYAAABdoNd6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtUVOe9PvBnuOMNhMFAig6o6AwagQYYTBPkeFKwi4Lx\nkqXYxq6YdI1pE4iS+8+zAsk6aS7NKZaVGpLW9ljpiV21iZfEEFzpdhITYDTGJChRE6hJDIoaEeTi\nAO/vj50ZGAfmAox7YJ7PWnvN7Mu75zsoz7y8e+/ZKiGEABERjWt+ShdARESex7AnIvIBDHsiIh/A\nsCci8gEMeyIiH8CwJyLyAU7D3mg0QqfTISEhAeXl5UNuZzKZEBAQgJ07d1qXxcXFYcGCBUhJSUF6\nevroVExERG4LcLZBUVERKioqoNFokJOTg4KCAqjVapttent78eijj2LJkiU2y1UqFSRJQkRExOhW\nTUREbnHYs29tbQUAZGZmQqPRIDs7G7W1tXbblZeXY+XKlYiKirJbx2u2iIiU57BnbzKZoNVqrfOJ\niYmoqalBbm6uddk333yDXbt24d1334XJZIJKpbKuU6lUWLx4MeLj47Fu3Trk5+fb7H/gtkRE5Dp3\nO9IjPkD74IMP4tlnn4VKpYIQwqaAgwcP4ujRo/jNb36DjRs3orm5edCCRzSNxj6cTE8++aTHX4N1\net80FurEL36heA3j5Wc5luocDodhn5aWhoaGBut8fX09MjIybLY5fPgwVq9ejfj4eOzcuRO/+tWv\nsHv3bgBATEwMAECn0yE/Px979uwZVpFERDQyDsM+LCwMgHxGTlNTE6qrq6HX6222+fLLL9HY2IjG\nxkasXLkSW7ZsQX5+Pjo6OtDW1gYAaGlpQVVVld0BXCIiuj6cno1TVlYGg8EAs9mMwsJCqNVqVFRU\nAAAMBsOQ7Zqbm7F8+XIAQGRkJIqLizF9+vRRKvv6ysrKUroEl7DO0TUm6kxOVroCl4yJnyXGTp3D\noRLDHQAajRf/fpx/hDsBlHsLRIpSSRLEOA4oGtxwspNX0BIR+QCGPRGRD2DYExH5AIY9EZEPYNgT\nEfkAhj0RkQ9g2BMR+QCGPRGRD2DYExH5AIY9EZEPYNgTEfkAxcO+oEDpCoiIxj/Fw/6115SugIho\n/FM87ImIyPMY9kREPsBp2BuNRuh0OiQkJKC8vHzI7UwmEwICArBz50632xIRkWc5DfuioiJUVFRg\n//79eOmll3D+/Hm7bXp7e/Hoo4/a3XbQlbZEROR5DsO+tbUVAJCZmQmNRoPs7GzU1tbabVdeXo6V\nK1ciKirK7bZEROR5Du9BazKZoNVqrfOJiYmoqalBbm6uddk333yDXbt24d1334XJZIJKpXK5rawE\nJSXys6ysrHF9D0giouGQJAmSJI1oH05vOO7Mgw8+iGeffdZ6T0T37ynbH/ZERGTv2o5waWmp2/tw\nGPZpaWl4+OGHrfP19fV24/KHDx/G6tWrAQDnz5/Hvn37EBgYiEWLFjltS0RE14fDMfuwsDAA8lk1\nTU1NqK6uhl6vt9nmyy+/RGNjIxobG7Fy5Ups2bIF+fn5LrUlIqLrw+kwTllZGQwGA8xmMwoLC6FW\nq1FRUQEAMBgMbrclIqLrTyXcH2QfvRdXqQAIjKgClQoj2wHR2KWSJAie1OBzLMdI3cEraImIfADD\nnojIBzDsiYh8gOJhz2O2RESep3jY9/QoXQER0fjHsCci8gEMeyIiH6B42Pf2Kl0BEdH4p3jYm828\nJoqIyNMUD3uVCujrU7oKIqLxTfGwDwyUe/dEROQ5iod9QAAP0hIReZriYc+ePRGR5zHsiYh8AMOe\niMgHMOyJiHyA07A3Go3Q6XRISEhAeXm53fpdu3YhKSkJycnJyM3Nhclksq6Li4vDggULkJKSgvT0\n9EH3z7AnIvI8p7clLCoqQkVFBTQaDXJyclBQUGBze8Hbb78dS5cuBQAcOHAAxcXFMBqNAOS7qUiS\nhIiIiCH3HxTEsCci8jSHPfvW1lYAQGZmJjQaDbKzs1FbW2uzzcSJE222DwkJsVnv7NZZgYHA1atu\n1UxERG5y2LM3mUzQarXW+cTERNTU1CA3N9dmu9dffx0bNmxAe3s7Dh06ZF2uUqmwePFixMfHY926\ndcjPz7d7jfPnS/CHPwA33ghkZWUhi/fTJCKyIUkSJEka0T6cDuO4YtmyZVi2bBl27NiBZcuW4ciR\nIwCAgwcPIiYmBsePH0deXh7S09MRHR1t01ajKcHatcCtt45GJURE48+1HeHS0lK39+FwGCctLQ0N\nDQ3W+fr6emRkZAy5/apVq3DmzBl0dnYCAGJiYgAAOp0O+fn52LNnj10bjtkTEXmew7APCwsDIJ+R\n09TUhOrqauj1epttvvjiC+u4/FtvvYWbb74ZoaGh6OjoQFtbGwCgpaUFVVVVWLJkid1rBAVxzJ6I\nyNOcDuOUlZXBYDDAbDajsLAQarUaFRUVAACDwYCdO3di27ZtCAwMREpKCp5//nkAQHNzM5YvXw4A\niIyMRHFxMaZPn263f4Y9EZHnqYSz02U8+eIqFZYvF1izBlixYtg74Rfik89SSRIET2rwOSqVyumZ\njtdS/Apa9uyJiDzPK8K+u1vpKoiIxjevCHv27ImIPEvxsA8OZtgTEXma4mHPYRwiIs9TPOzZsyci\n8jyvCHv27ImIPEvxsOcwDhGR5yke9iEhHMYhIvI0xcM+OBjo6lK6CiKi8U3xsA8JYdgTEXma4mHP\nA7RERJ6neNizZ09E5HmKhz3H7ImIPE/xsA8J4TAOEZGneUXYf38XQyIi8hCnYW80GqHT6ZCQkIDy\n8nK79bt27UJSUhKSk5ORm5sLk8nkcluAPXsioutCOJGcnCwOHDggmpqaxNy5c0VLS4vN+vb2dutz\nSZLEbbfd5nJbAOKTT4SYN89ZFQ44fwtE4xb+9S+lSyAFuBDddhz27FtbWwEAmZmZ0Gg0yM7ORm1t\nrc02EydOtNk+JCTE5bYAz8YhIroeHN5w3GQyQavVWucTExNRU1OD3Nxcm+1ef/11bNiwAe3t7Th8\n+LBbbbdsKcG5c0BJCZCVlYUs3k+TiMiGJEmQJGlE+3AY9q5atmwZli1bhh07duCOO+7AkSNHXG77\nxBMl+N//lcOeiIjsXdsRLi0tdXsfDodx0tLS0NDQYJ2vr69HRkbGkNuvWrUKZ86cQWdnJ1JTU11q\ny2EcIiLPcxj2YWFhAOSzapqamlBdXQ29Xm+zzRdffAH5eAHw1ltv4eabb0ZoaCjCw8OdtgX6w/77\nXRARkQc4HcYpKyuDwWCA2WxGYWEh1Go1KioqAAAGgwE7d+7Etm3bEBgYiJSUFDz//PMO29oVEAD4\n+QFms/zd9kRENPpUQijXp1apVBBCYPJk4Ouvge//kHB3J/yzgHyWSpIgeFKDz7FkpzsUv4IWAEJD\nOW5PRORJXhP2/MoEIiLPYdgTEfkAhj0RkQ/wirDnufZERJ7lFWHPnj0RkWcx7ImIfADDnojIBzDs\niYh8gFeEPW9NSETkWV4R9uzZExF5lteEPU+9JCLyHK8Je/bsiYg8h2FPROQDGPZERD6AYU9E5AOc\nhr3RaIROp0NCQgLKy8vt1ldWViIpKQlJSUlYs2YNTpw4YV0XFxeHBQsWICUlBenp6UO+BsOeiMiz\nnIZ9UVERKioqsH//frz00ks4f/68zfqZM2fCaDTi6NGjyMnJwdNPP21dp1KpIEkSjhw5grq6uiFf\ng2FPRORZDsO+tbUVAJCZmQmNRoPs7GzU1tbabLNw4ULrjclzc3Nx4MABm/Wu3DqLYU9E5FkObzhu\nMpmg1Wqt84mJiaipqUFubu6g27/yyivIy8uzzqtUKixevBjx8fFYt24d8vPz7dqUlJSgqQk4fhyQ\npCxk8X6aREQ2JEmCJEkj2ofDsHfH/v37sX37dnzwwQfWZQcPHkRMTAyOHz+OvLw8pKenIzo62qZd\nSUkJamrksGfOExHZy8qy7QiXlpa6vQ+HwzhpaWloaGiwztfX1yMjI8Nuu08++QTr16/H7t27ER4e\nbl0eExMDANDpdMjPz8eePXsGfZ0JE4CODrdrJyIiFzkMe8tYvNFoRFNTE6qrq6HX6222OX36NFas\nWIHKykrMnj3buryjowNtbW0AgJaWFlRVVWHJkiWDvg7DnojIs5wO45SVlcFgMMBsNqOwsBBqtRoV\nFRUAAIPBgKeeegoXL17E+vXrAQCBgYGoq6tDc3Mzli9fDgCIjIxEcXExpk+fPuhrTJwIXLkyWm+J\niIiupRKunC7jqRdXqSCEwOXLQGwscPnysHYCKPcWiBSlkiQIHuzyOZbsdIdXXEE7YYLcs2dmExF5\nhleEfUCAPHV3K10JEdH45BVhD8jj9jxIS0TkGV4T9pahHCIiGn1eE/bs2RMReY5XhT179kREnuE1\nYc8Lq4iIPMdrwp49eyIiz/GasOcBWiIiz/GasJ80iWFPROQpXhX27e1KV0FEND4x7ImIfADDnojI\nB3hN2PNsHCIiz/GasGfPnojIcxj2REQ+wGnYG41G6HQ6JCQkoLy83G59ZWUlkpKSkJSUhDVr1uDE\niRMutx2IYU9E5DlOw76oqAgVFRXYv38/XnrpJZw/f95m/cyZM2E0GnH06FHk5OTg6aefdrntQAx7\nIiLPcRj2ra2tAIDMzExoNBpkZ2ejtrbWZpuFCxdab0yem5uLAwcOuNx2IIY9EZHnOLzhuMlkglar\ntc4nJiaipqYGubm5g27/yiuvIC8vz622JSUlAIBz54CzZ7MAZLn/LoiIxjFJkiBJ0oj24TDs3bF/\n/35s374dH3zwgVvtLGHf2Ai89dZoVUNENH5kZWUha8CN5UtLS93eh8NhnLS0NDQ0NFjn6+vrkZGR\nYbfdJ598gvXr12P37t0IDw93q60Fh3GIiDzHYdhbxuKNRiOamppQXV0NvV5vs83p06exYsUKVFZW\nYvbs2W61HYhhT0TkOU6HccrKymAwGGA2m1FYWAi1Wo2KigoAgMFgwFNPPYWLFy9i/fr1AIDAwEDU\n1dUN2XYoISFATw9w9SoQFDQab42IiCxUQgih2IurVBj48tHRwEcfATfe6NZOAOXeApGiVJIEMWAs\nl3zDtdnpCq+5ghYApk2Tz8ohIqLRxbAnIvIBXhX2N9zAsCci8gSvCnv27ImIPINhT0TkAxj2REQ+\ngGFPROQDGPZERD7Aq8I+Ohr49lulqyAiGn+8KuxjY4Hz54HOTqUrISIaX7wq7P39gfh44NQppSsh\nIhpfvCrsAWDOHGDAbWyJiGgUMOyJiHwAw56IyAd4XdgnJAAnTypdBRHR+OI07I1GI3Q6HRISElBe\nXm63vqGhAQsXLkRISAhefPFFm3VxcXFYsGABUlJSkJ6e7lJBc+YAn3/uYvVEROQSp3eqKioqQkVF\nBTQaDXJyclBQUGBzx6nIyEiUl5fjjTfesGurUqkgSRIiIiJcLig6Wr5j1blz8kVWREQ0cg579q2t\nrQCAzMxMaDQaZGdno7a21mabqKgopKamIjAwcNB9uHs3FZUKmD8fqK93qxkRETngsGdvMpmg1Wqt\n84mJiaipqUFubq5LO1epVFi8eDHi4+Oxbt065Ofn221TUlJifZ6VlYWsrCzMnw989hnwH//h4rsg\nIhrHJEmCJEkj2ofTYZyROHjwIGJiYnD8+HHk5eUhPT0d0dHRNtsMDHuL+fOBo0c9WRkR0dhh6Qhb\nlJaWur0Ph8M4aWlpaGhosM7X19cjIyPD5Z3HxMQAAHQ6HfLz87Fnzx6X2t10k9yzJyKi0eEw7MPC\nwgDIZ+Q0NTWhuroaer1+0G2vHZvv6OhAW1sbAKClpQVVVVVYsmSJS0XNmyeHvZvD/URENASnwzhl\nZWUwGAwwm80oLCyEWq1GRUUFAMBgMKC5uRlpaWm4fPky/Pz8sHnzZhw7dgznzp3D8uXLAchn7BQX\nF2P69OkuFRUZCUyaBHz1FTBjxgjeHRERAQBUwt3TZUbzxVWqIc/WyckBfvELYM0apzvhnwDks1SS\nBDFgLJd8g6PsHIrXXUFr8cgjwGOPAZcvK10JEdHY57Vh/5//KffuH3tM6UqIiMY+rw17AHjhBWDH\nDuD0aaUrISIa27w67MPDgfx84PXXla6EiGhs8+qwB4Dly4F//lPpKoiIxjavD/sf/1i+mvbsWaUr\nISIau7w+7ENCgCVLgN27la6EiGjs8vqwBziUQ0Q0UmMi7H/yE+DgQeDSJaUrISIam8ZE2E+eDGRk\nAO+9p3QlRERj05gIewBYuBD44AOlqyAiGpvGTNjfcgvw4YdKV0FENDaNmbDX64FDhwCzWelKiIjG\nnjET9uHhQFwc8MknSldCRDT2jJmwBzhuT0Q0XGMq7DluT0Q0PE7D3mg0QqfTISEhAeXl5XbrGxoa\nsHDhQoSEhODFF190q627brmFPXsiouFweqeqlJQUbN68GRqNBjk5OXj//fehVqut61taWvDvf/8b\nb7zxBqZOnYri4mKX27p7txUhALUa+PRT4MYbrTvhnarIZ/FOVb5p1O9U1draCgDIzMyERqNBdnY2\namtrbbaJiopCamoqAgMD3W7rLpVKHrfnUA4RkXsc3nDcZDJBq9Va5xMTE1FTU4Pc3FynO3a1bUlJ\nifV5VlYWspz0Uixhv2KF0xKIiMYFSZIgSdKI9uEw7K+HgWHviltukW9VKITc03fFpUvAX/8KpKbK\n5+v7janD0kTk667tCJeWlrq9D4exl5aWhoaGBut8fX09MjIyXNrxSNo6kpEBdHQAhYVAb6/jbYWQ\nQz4xEThwAPjlL4EZM4CiIvl7dvr6RlwOEdGY4DDsw8LCAMhn1TQ1NaG6uhp6vX7Qba89WOBOW3eE\nhspBfeyY46GcTz8FFi0CNm8G3ngD+Mc/gM8+A6qr5YO8v/41EBsL3H8/jwEQkQ8QTkiSJLRarZg1\na5bYvHmzEEKIl19+Wbz88stCCCG+/fZbERsbK6ZMmSLCw8PF9OnTRVtb25BtB3Lh5YfU3S3E2rVC\nCED89a9C/OMfQuzdK8T+/UJs3ChEVJQQW7YI0dMz9D4aGoT47/8WYvp0IX71KyGuXBl2OUSKwL/+\npXQJpIDhZKfTUy89aTinDw0kBKDyU+HnPxPo7AS6uoDOTmDuXOCpp4CoKNf2c+mS3MM/dEge9klL\nG3ZJRNcVT730TcPJzjEd9t/vZNTOs9+xA3jgAXl6/HEgQPHD10SOMex906ifZ+9rVq0CPvpIPpib\nmQl88YXSFRERjQ6G/TViY4F33pGDPyMDePVVXqBLRGMfw34Qfn7y6ZmSBPzhD8AddwAtLUpXRUQ0\nfAx7B+bNA2pr5QO+qanyEA8R0VjEsHciKAh4/nngxReBnBzgb39TuiIiIvfxfBMXrVwJzJkjD+kc\nOQI8+yzg7690VUREruGpl266cEE+eBsQAPzf/wFTp163l7bT2QnU1wNHj8rT55/LP4qgICAwsP9x\n4PPBljlbHxkJ6HRARIRy75UGx1MvfdNwspM9ezdFRgJvvw088giQni5/FcO8eZ59TSGAM2f6Q90y\nNTXJf20sWAAkJcnDTP7+8k3Zr16VHy3TwPmBzzs6Bl8+cP7sWeD4cWDCBDn0ExNtH6OjXf9SOiJS\nBnv2I7BtG/DQQ/LpmUuXjs4+u7vl7/05elS+ubol2FUqOdAHTjqd3Pu+HiwfOMeOycE/8LGnZ/AP\ngRkz+A2jnsaevW/iFbQKMJmA5cuBe+8F/uu/XA+3ri7g5El56OXzz+XgPHoUOHUKmDnTNtQXLABi\nYry399zSItd/7YfApUvymUzXfgjMmsWrk0cLw943MewV0twsfwPntGlyb3/yZHm5pTdsCXTL1NAA\nfPstEB8vh+HcuYBWK4f6vHlASIiib2fUtLbK7/XaD4EzZ+TAv/ZDYM6c8fPerxeGvW9i2Cuou1v+\nTp2DB4Hk5P5gnzChP9AtoT53rhz0vtq77ezs/2tm4IfAl18CcXHyNQ2pqfIX0qWkyD9DGhzD3jcx\n7BUmBPD660B7e3+4h4crXdXYYTbLwX/okDw8duiQfLZRQkJ/+KelATfddP2OVXg7hr1vYtjTuNPd\nLR+otoS/ySR/Qd38+f3hn5Ym/8Xki9c9MOx9E8OefEJ7O/Dxx3LwW6bmZnnIxxL+qanycQFvPag9\nWhj2vskjYW80GmEwGNDT04PCwkI88MADdts8/vjj2LFjB6ZOnYrKykpotVoAQFxcHKZMmQJ/f38E\nBgairq5uxAXbvwOGPQHffQccPtwf/ocOyR8KluEfy+MPfjC+PgAY9r7JI2GfkpKCzZs3Q6PRICcn\nB++//z7UarV1fV1dHTZu3Ijdu3ejqqoKlZWV2Lt3LwAgPj4ehw8fRsQQl14y7MmTmpttx/9NJnmo\nZ2D4p6XJ9yQeqxj2vmnUr6BtbW0FAGRmZgIAsrOzUVtbi9zcXOs2tbW1WLlyJSIiIlBQUIBNmzbZ\n7EPBUSLycdHRwE9/Kk+A3Cc4fbo/+H/7W/mvgalTbcP/hz8EwsKUrZ1otDkMe5PJZB2SAYDExETU\n1NTYhH1dXR3uuusu63xUVBS+/PJLzJw5EyqVCosXL0Z8fDzWrVuH/Px8u9coKSmxPs/KykIWeynk\nISoVoNHI04oV8rK+PvlCNsvwz6ZN8sVtN94on02VkCCf/5+QIE+xsbwqmK4/SZIgSdKI9jHiM72F\nEEP23g8ePIiYmBgcP34ceXl5SE9PR3R0tM02A8Oe6Hrz85PDfM4c4Gc/k5f19MjXAZw8CZw4Iff+\nX3tNnv/uO/kK54EfAJbJm69yprHt2o5waWmp2/twGPZpaWl4+OGHrfP19fVYsmSJzTZ6vR7Hjh1D\nTk4OAKClpQUzZ84EAMTExAAAdDod8vPzsWfPHvzyl790u0ii6ykgQL6SebAvuGtvl/8SOHlSng4e\nBP7yF/lDoaPD/gPAMkVF8YOAlOUw7MO+H7g0Go2YMWMGqqur8eSTT9pso9frsXHjRqxduxZVVVXQ\n6XQAgI6ODvT29mLy5MloaWlBVVUVNmzY4KG3QXR9TJokXyGdnGy/rrW1/0Pg5Elg/35gyxb5eV+f\n/QeA5a8DJb8mm3yH02GcsrIyGAwGmM1mFBYWQq1Wo6KiAgBgMBiQnp6OW2+9FampqYiIiMD27dsB\nAM3NzVi+fDkAIDIyEsXFxZg+fboH3wqRssLC+r/q4VoXLvR/CJw4Aezd2z8fFDT4h0BCQv/3LBGN\nFC+qIlKQEMC5c/IHwMC/Ck6elIeLJk+2/wBISABmz5a/M4inXvomXkFLNI709cnfEDrwA8DyodDY\nKN9I55vtEu7elgV/f3l7IQZ/dLTO2aNSbQf+WqtU/cc8Bj46WzacNt68n9JSYPFihj2Rz+jtBb76\nCohvkvDqqSwIIf8q+Pn1Pw587u7jSNqO5usD8q+35Vd84KOzZcNp4+370WotB/sZ9kQ+hcM4vmk4\n2cnLQ4iIfADDnojIBzDsiYh8AMOeiMgHMOyJiHyAj97ymohodAgh0Cf60NPXg56+HvSK3v7nfb2D\nLne0zlGbRZpFiJ8aP6w6GfZEPsoSUn2iD72iV37s6/Wa+d6+XmvYDRmMwwhQV8PVnTZ+Kj/4q/wR\n4Bdgnfz9+ucHrhu43NG6wfanU+sY9kTXkxACPX09uNp7FVd7r6K7t7v/eU/3kMtc2tZJu4HLoC3D\nrN/PcjlIBz4XEFBBBX8/f2tY+an8vGreYWiq/BHsH+xWuHoikC31ejteVEVep7evF+Y+s02ouTIN\nDMnR2N5Z0Pqp/BDkH4TggGD50V9+HLhs4PLBltks93N9X5blKQ2XcHJB7LCDVcXvXR6TRv22hOQb\nhBAw95nR1dOFrp4udPd09z/v7bYuu/b5tY+W7Qdb7+q6q71X0Sf6rIE2WOC5Mg0MRcs0IXCCy9sH\n+gci2D94yKAN8g+Cv5+/0v90QIOE2RGzla6CxgCGvcLMvWZcMV/BlatXhgxDZyE7WHheG9ZDhbhl\neYBfAIIDghESEIKQgBAE+w94HhBsDT7L8muXWR6nBE8ZdLll387WBfkHIcAvgD1OolHGsHeBEAKd\nPZ1ov9qO9qvtuHL1Sv9z8xW75ZZllhAf+GjZzrJMQGBi4ERMCJzgMEQdPYYGhCI8ONxhWA+2fODr\neUUvlYg8ZlyGvSWc27rb0H61HW1X29DW3Ya2q2243H0Zbd3y4+Wr/c8tgd12ta3/+fft20+0I3h2\nMCYFTcKkoEmYGDix/3nQRJvlE4MmIjI0EjPCZtgss7SxPLe0C/IPGrX3LUnSmLhhO+scRR9/DHh7\njRgjP0uMnTqHw2nYG41GGAwG9PT0oLCwEA888IDdNo8//jh27NiBqVOnorKyElqt1uW2rjD3mnHu\nyjmcvXIWZ9vP4uyVszjfcR4tHS14DsDS15bifMd563Sp6xKC/IMwOWgyJgVNwuRg+XFK8BRMCZ6C\nyUGTrc9/MPkH0Kq1dttOCppkXfY/v/kflP4/92/we72Nlf+orHMUffyx0hW4ZEz8LDF26hwOp2Ff\nVFSEiooKaDQa5OTkoKCgAGq12rq+rq4O7733Hg4dOoSqqio89NBD2Lt3r0ttB7rYeRGfnfsMx1qO\n4fMLn+PUxVP4qvUrnGk7g++6voN6gho3TLwBN0y6AdMmTkPUhChETYgCANydfDfUE9RQT1AjMjQS\n4SHhCPQPHI2fDwBw/JiIxjyHYd/a2goAyMzMBABkZ2ejtrYWubm51m1qa2uxcuVKREREoKCgAJs2\nbXK5LQA8tv8x/PP4P9Hc3oz50+YjMSoRcyLnYJFmETRhGsRMjkHUhCgHY8pP4A7tHe6/cyIiXyIc\nqK6uFqtXr7bOb9myRWzatMlmm5///OeiqqrKOq/X68WpU6dcaguAEydOnDgNY3LXiA/QCiHsTu53\nddhDweu5iIh8isNrfNPS0tDQ0GCdr6+vR0ZGhs02er0ex44ds863tLRg5syZSE1NddqWiIiuD4dh\nHxYWBkA+q6apqQnV1dXQ6/U22+j1euzcuRMXLlzA3/72N+h0OgBAeHi407ZERHR9OB3GKSsrg8Fg\ngNlsRmGFke+yAAAF+0lEQVRhIdRqNSoqKgAABoMB6enpuPXWW5GamoqIiAhs377dYVsiIlKA26P8\no+TAgQNCq9WK2bNni9///vdKleHQ6dOnRVZWlkhMTBSLFi0SlZWVSpfkUE9Pj0hOThY//elPlS5l\nSO3t7WLt2rUiISFB6HQ68eGHHypd0qBeeeUVsXDhQvHDH/5QFBUVKV2O1d133y2mTZsm5s+fb112\n+fJlkZ+fL6ZPny6WLl0q2traFKxw8BofeughodVqRUpKiigqKhIdHR0KVigbrE6L3/72t0KlUokL\nFy4oUJmtoercunWr0Gq1IjExUTzyyCNO96NY2CcnJ4sDBw6IpqYmMXfuXNHS0qJUKUP69ttvxZEj\nR4QQQrS0tIj4+Hhx+fJlhasa2osvvijWrFkj8vLylC5lSMXFxWLTpk2is7NTmM1mcenSJaVLsnPh\nwgURFxcn2tvbRW9vr/jJT34i3n77baXLEkIIYTQaxUcffWTzi//cc8+J+++/X3R1dYlf//rX4oUX\nXlCwwsFrfOedd0Rvb6/o7e0V9957r/jjH/+oYIWyweoUQu7k5eTkiLi4OK8I+8Hq/PTTT0VGRoY4\nceKEEEKIc+fOOd2PIl/CPPAcfI1GYz0H39tER0cjOTkZAKBWqzFv3jwcOnRI4aoG9/XXX+Ott97C\nvffe69VnOe3fvx9PPPEEQkJCEBAQYD0u5E1CQ0MhhEBrays6OzvR0dGBqVOnKl0WAOC2226zq6Wu\nrg733HMPgoODsW7dOsV/lwar8cc//jH8/Pzg5+eHnJwcHDhwQKHq+g1WJwBs3LgRzz//vAIVDW6w\nOvft24d77rkHCQkJAICoqCin+1Ek7E0mk/UrFQAgMTERNTU1SpTislOnTqG+vh7p6elKlzKoDRs2\n4IUXXoCfn/feROHrr79GV1cX7rvvPuj1ejz33HPo6upSuiw7oaGh2LJlC+Li4hAdHY0f/ehHXvvv\nDtj+Pmm1WtTV1SlckWOvvvoq8vLylC5jULt27UJsbCwWLFigdCkOvfPOO/jss8+QmpqKe++91+aM\nyKF4bzJ4kba2NqxatQq/+93vMHHiRKXLsbN3715MmzYNKSkpXt2r7+rqwokTJ7BixQpIkoT6+nr8\n/e9/V7osOy0tLbjvvvtw7NgxNDU14cMPP8Sbb76pdFlD8uZ/82s99dRTmDx5Mu68806lS7HT0dGB\nZ555BqWl/d+D5a0/266uLly8eBHvvfceli5divvvv99pG0XC3pXz972F2WzGihUrcNddd2Hp0qVK\nlzOoDz74ALt370Z8fDwKCgrw7rvvYu3atUqXZWf27NmYO3cu8vLyEBoaioKCAuzbt0/psuzU1dUh\nIyMDs2fPRmRkJO68804YjUalyxpSWloajh8/DgA4fvw40tLSFK5ocH/5y19QVVVlc8aeN/niiy/Q\n1NSEpKQkxMfH4+uvv8bNN9+Mc+fOKV2anYyMDKxatQqhoaHIy8tDQ0OD07+SFQl7V87f9wZCCNxz\nzz2YP38+HnzwQaXLGdIzzzyDr776Co2NjXjttdewePFibNu2TemyBpWQkIDa2lr09fXhzTffxO23\n3650SXZuu+02HDp0CBcvXkR3dzf27duH7Oxspcsakl6vx9atW9HZ2YmtW7d6Zcfp7bffxgsvvIDd\nu3cjJCRE6XIGddNNN+Hs2bNobGxEY2MjYmNj8dFHH2HatGlKl2Zn4cKF2LdvH4QQqK2txaxZs5z/\nXEf/2LFrJEkSWq1WzJo1S2zevFmpMhx67733hEqlEklJSSI5OVkkJyeLffv2KV2WQ5IkefXZOJ9/\n/rnQ6/UiKSlJFBcXi/b2dqVLGtSf//xnkZmZKVJTU8WmTZtEb2+v0iUJIYRYvXq1iImJEUFBQSI2\nNlZs3brV6069tNQYGBgoYmNjxZ/+9Ccxe/ZsMWPGDOvv0X333adojQPrHPizHCg+Pt4rzsYZrM6e\nnh5hMBiEVqsVd9xxh6irq3O6H0VvOE5ERNcHD9ASEfkAhj0RkQ9g2BMR+QCGPRGRD2DYExH5AIY9\nEZEP+P+SSFsWJ0cZLwAAAABJRU5ErkJggg==\n"
      }
     ],
     "prompt_number": 39
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "bins, CDF = fit2.cdf()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 51
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "len(bins)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 52,
       "text": [
        "46"
       ]
      }
     ],
     "prompt_number": 52
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fit2 = powerlaw.Fit(data, xmin =.794)\n",
      "fit2.alpha\n",
      "fit2.sigma\n",
      "fit2.plot_ccdf()\n",
      "fit2.D"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 43,
       "text": [
        "0.092103799365054417"
       ]
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEHCAYAAABGNUbLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHuhJREFUeJzt3Xt8VOWdx/FPuCkaTaVgREFQjJJIuKybBtGEQWjCiimK\noI0aNCESQRJQFxFomwBqC7UKBgwgyKp9FV3XCMjVRJoEqEQ0YFOURYTZrVjXAiWlyE1z9o9HUiIB\nZiYzc2bmfN+vV17NHJI5v6cTf/PM77lFWZZlISIijtLC7gBERCT4lPxFRBxIyV9ExIGU/EVEHEjJ\nX0TEgZT8RUQcSMlfRMSBlPxFRByoVaCe+NixY0yePJkjR44wdOhQBg8eHKhbiYiIlwLW89+0aRNJ\nSUmUlJRQWloaqNuIiIgPvEr+OTk5xMbGkpiY2Oh6VVUV8fHxxMXFUVxcDEBtbS3dunUD4MiRI34K\nV0RE/MGr5J+dnc3atWtPuz5+/HgWLFhAeXk58+bNY9++ffTs2ZPdu3cDcMEFF/gnWhER8Quvav4p\nKSm43e5G1+rq6gBITU0FIC0tjerqagYNGsTUqVPZtGkTw4YNa/L5oqKifAhZRESauydns2v+W7Zs\noXv37g2PExIS2Lx5M+eddx7PPPMMxcXFpKenn/H3Lctq1ldhYWGzf66pf/Pk2qmPc3MLmT7dIiXF\nok2bQgYMsBg4sJDa2vBs39naevJ7T2MLhfZ50p5z/X8QqPZ527ZQaV+gXjt/tC+c/jZ9aZ8/tCwq\nKiry5hcOHjzI0qVLGTt2LAC7d+9m69atDB8+HIAPPviAQ4cOccstt5zzuaZNm9bwfdeuXb0JoxFP\nf/dsP9fUv3ly7eTjiy+G++/vSnY23HAD9O3blfLyCubNc/Huu9C1K1x5pUdhehW3pz/nS/vO1NaT\n31dUVOByuTyK7WyC1b5ztef73wezfd62ranrdrQvUK9dU9e9aV+4/W02da2p9m3bto21a9dSWVmJ\nl6n7dJaX9uzZY/Xo0aPh8cGDB63evXs3PB43bpy1cuVKj57Lh9uHlcLCQuvIEcsqKbGsrl0ty+Wy\nrLIyy6qvtzsy/ygsLLQ7hIBS+8JXJLfNsvyTO5td9omJiQHMjB+3201ZWRnJycke/35RUREVFRXN\nDSMkuVwuzj8fHnoIdu6E7GzIz4e+feG11+DoUbsjbB5/9KxCmdoXviK1bRUVFc3v8X8n6rt3EY9k\nZmZSWVnJ/v37ufTSS5k+fTrZ2dlUVlby0EMPceLECQoKCigoKPDs5lFRfqtfhYtvv4Vly2D+fNi6\nFe69F3Jz4XuzZ0VEzsgfudOr5O9vTkz+p9qzB5YsMV+XX27eBDIzITra7shEJJT5I3favrdPJJd9\nzuWqq2D6dHC7oagI1qwxg8IPPwy1tXZHJyKhxrayj785vefflL17YdEiePFFM0PooYdg+HA4/3y7\nIxORUKGyTwT75htYuRJKSqCmBh54APLy4Jpr7I5MROymsk8Ea9UKbr8d1q2D994z1268EdLT4a23\nzJuDiDiLyj4OdfQo/Nd/mU8D//M/MHo0PPggdOxod2QiEkwR0fMXz51/Ptx3H2zaBKtWwRdfQEIC\n3H03VFaC3kdFxFO2J3+VfXzTq5dZK+B2Q0oKjBkDPXqYaaP19XZHJyKBoLKPnMayoKICJk+GNm1g\n4UI4Zb89EYkgKvtIg6goGDDAlITuugtuvhmmTYNjx+yOTERCkZJ/hGnZEsaNM1tHfPgh9OkDGzfa\nHZWIhBqVfSKYZUFpKRQUQP/+5pNAXJzdUYlIc0VE2UcDvoETFQV33gn//d9w/fVmncCDD8Kf/2x3\nZCLiCw34ik8OHIBnnjGzhEaONIPDsbF2RyUi3oqInr8ET7t28PTT8PHHZjpoQgL87Gdw8KDdkYlI\nsCn5O9Bll8Hzz5s9g/buhWuvhVmz4Ouv7Y5MRIJFyd/BunQxi8IqK+H9981gcEkJnDhhd2QiEmi2\nJ38N+NovPt7sGbR8udk0LjHRnC0gIqFFA74SMJYFq1fDI4+YTwLPPgvXXWd3VCJyKg34it9FRcGQ\nIfCnP8Ett8BNN8Fjj2lQWCTSKPlLk9q0MUl/+3b4+9/NPkHPPQdHjtgdmYj4g5K/nFVsrDlS8p13\noKrKnCT2wgvaM0gk3Cn5i0d69jSDwcuXm+Mlr7sOFi/WzCCRcKUBX/HJH/4AU6eaTwBvvw0//KHd\nEYk4hwZ8xTb9+sH69WbDuJtvNsdKikj4aGV3AEVFRbhcLlwul92hiJeiouCXvzRnCN90k5ki2rOn\n3VGJRK6Kigq/rYtS2Uf84vXXIT8f/vM/Qe/jIoGlso+EjLvvhtdeM6eIvfGG3dGIyLmo5y9+tW0b\n3HYb3HEHPPUUXHyx3RGJRB71/CXk9O4Nf/wjHD1qtox+802zZYSIhBb1/CVgNmyAvDzo1g3mzjW7\niIpI86nnLyEtJcWUgZKT4YYb4De/gW++sTsqEQH1/CVIPv3UfAo4dAgWLYJeveyOSCR8qecvYSMu\nDt59F8aMgR//GKZM0SZxInYKaPLfs2cPubm5jBgxIpC3kTARFQU5OWZAeNcu0/uvrLQ7KhFnCmjy\nv+qqq1i0aFEgbyFh6LLLzGKwX/8a7rsPRo/WeQEiweZR8s/JySE2NpbExMRG16uqqoiPjycuLo7i\n4uKABCiRa+hQc2hMq1Zw/fVQWmp3RCLO4VHyz87OZu3ataddHz9+PAsWLKC8vJx58+axb98+Xn31\nVR555BG++OILvwcrkScmxpwP8PrrZpfQYcNAfzoigefRxm4pKSm43e5G1+rq6gBITU0FIC0tjerq\narKyssjKygLgwIEDTJkyhW3btjFz5kwmTZp02nOfehixNnhzrptvhq1b4emnzVjAk0/Cgw9CC01J\nEPHrhm4neTzV0+12k5GRQW1tLQDl5eUsXryYpUuXAjB//nz27t3LjBkzPL+5pnpKE/70J8jNhdat\nYeFCiI+3OyKR0BIRUz2Lior8/o4m4a1HD9i0CX76U0hNhWnTdGykCJhPAKdWS5rD555/XV0dLpeL\nrVu3ApCfn8/gwYMZMmSI5zdXz1/O4fPP4eGHzSKxhQtNeUjE6Wzt+cfExABmxo/b7aasrIzk5ORm\nBSPyfZ06wbJlMGOG2Tb68cfh22/tjkok/HmU/DMzM+nXrx87d+6kc+fOLFmyBIDZs2eTl5fHoEGD\nGDt2LO3bt/c6AJV95FyiouDOO6G2FmpqICMDvptvIOIotpR9AkFlH/HWiRPwyCNmq4gVK8y2ESJO\nowFfcZzWrc320OPHm/r/u+/aHZFI8KjnLwL8/veQmQk/+5kZFI6KsjsikeDwR+5U8pew9tlnMHw4\nXHIJFBebbSJEIp3KPuJ43brBli1mWwiXCx57DP7+d7ujEgkMlX1EmvDVV/DEE7BuHcyaBffco1KQ\nRCaVfUSa8N57ZgygXTtYsgQ6d7Y7IhH/ioiyj4i/3XijKQUNHGjODn7tNbsjEgk9tid/1fwlEFq2\nhMmTYc0aKCqCe+/VgTES/lTzF/HC11/DxImwciW8/LIZGBYJZ6r5i3hh9WqzVfTQofDUU2ZMQCQc\nqeYv4oVbb4Xt201JKD4eFi2C+nq7oxKxh+3JXzV/CaZLLjHbQ6xdCy+9BP36wYcf2h2ViGdU8xfx\ng/p6MwYweTLcdptZIKZTwyQcqOwj0gwtWkB2NnzyCVxxBQwYAIMHmxlCKgdJpFPPX+Q7R4/C66/D\n7Nlw5Ajk50NODrRta3dkIo1pto9IAFgWbNgAM2fCX/5iThK78kq7oxL5J5V9RAIgKsocHL9ypVkc\nlpxs3gxEIontyV+zfSRURUWZQeD/+A+zbfSCBXZHJE6n2T4iQfbpp2ZxWP/+MGcOtGljd0TiZCr7\niARJXBxs3gx795pZQR99ZHdEIs2j5C/ioYsvNoO/mZmQlgYPPACff253VCK+UfIX8UKLFjBuHOzc\nadYG9OoFU6ZAXZ3dkYl4R8lfxAcxMWZzuI8+MtNBr70WnnvOrA8QCQdK/iLN0KmTOS2srMxMB736\navMm8PXXdkcmcnZK/iJ+0LMnlJaarSE2bDAHy+tNQEKZ7clf8/wlkvTu3fhN4JprzDGSmtEs/qB5\n/iJhoroaRo0ynwRKSuDyy+2OSCKB5vmLhLjkZHNeQK9e5lPBSy/pU4CEBvX8RYLko4/MFtLt28PC\nhdC1q90RSbhSz18kjPTqZcpAAwbADTfAz38O//iH3VGJUyn5iwRR69bm5LBt22DPHrjuOjNVVIfH\nSLCp7CNio+pqePRRszjs2WfB5bI7IgkHOsxFJAJYFrzxBjz+OAwbBr/+NbRsaXdUEsqU/EUiyIED\nJvlfcgn89rdw4YV2RyShKuQHfJcvX87o0aPJycnh/fffD+StRMJeu3bwzjtw0UWm/PPll3ZHJJEs\nKD3/r776isLCQkpKShrfXD1/kdNYFsyYYdYErFoF119vd0QSaoLW88/JySE2NpbExMRG16uqqoiP\njycuLo7i4uIz/v7MmTPJy8trVqAiThEVBb/4BTz5pJkW+s47dkckkcijnv+GDRuIjo5m5MiR1NbW\nNlzv06cPc+bMoUuXLqSnp7Nx40bWrFlDTU0NEydOpGPHjkyaNIn09HQGDhx4+s3V8xc5q6oquOce\ns1L4ySchPt7uiCQU+CN3tvLkh1JSUnC73Y2u1X13ekVqaioAaWlpVFdXk5WVRVZWFgDPP/8869ev\n59ChQ+zatavJ3v+pmxS5XC5cmusm0iA11ZwfPHeuOT/41luhqEirg52moqLC7xtgelzzd7vdZGRk\nNPT8y8vLWbx4MUuXLgVg/vz57N27lxkzZnh+c/X8RTxWVwfPPAMvvAD33gtTp0JsrN1RiR1CfraP\nJ7Sls4hnYmLMQPAnn5jjJBMSzNqAv/7V7sgkWPy5pbPPyT8pKYkdO3Y0PN6+fTt9+/b1+nmKiopU\n6hHxwqWXwuzZZqO4w4ehe3fzKeDAAbsjk0BzuVz2J/+YmBjAzPhxu92UlZWRnJzs9fOo5y/im06d\nYN48qKkxvf9rr4XCQjh0yO7IJFCCfphLZmYmlZWV7N+/n0svvZTp06eTnZ1NZWUlDz30ECdOnKCg\noICCggLvbq6av4jf7N5tBoPLysybQG4utPJoSoeEG23vICKnqamBf/93s0J41iwYMsSsHZDIoQFf\nETnNv/wLvPuu2SBu0iQYONCcJ6x+VvjTGb4i4pFvvjHbRMyeDceOmSmi991nxgckfKnsIyIesSzY\nuhVefRWWLoUuXSArC/LyzAEzEl5U9hERj0RFmXLQc8/B55/DtGnw5pvmU8C339odnXhKZR8Rabaj\nRyEjA664wpSGWtjeFRRPRUTPX0Tscf75sGwZfPYZ5OdrQNhpbE/+KvuI2OfCC2HlSnj/fTMzSG8A\noU1lHxHxq/37zelhI0aYswQktAVtS2cRiWw//KFZGdy/P7RtCxMn2h2RBJqSv4gAcNllUF4Ot9wC\nx4+bzeIkcin5i0iDzp3N6WGDBpkdQ596SltDRCoN+IpIIx07QkUFrF0LEyZoEDiUaMBXRALu4EFz\nbOT118P8+dCypd0RyUma5y8iAfODH8A778CuXTByJJw4YXdE4k9K/iJyRtHRsHq1qf/37m3WBOjD\nemRQ2UdEzsmyYNUqc2ZwbKzZLvpf/9XuqJwrIso+GvAVCX1RUXDbbfDHP8I998BPfmL+d88euyNz\nFg34ioit/vEPePZZc07Aj38M48bBzTdrWmiwaD9/EbFVXR288grMnWs2inv4YXNgzIUX2h1ZZFPy\nF5GQUF9vjo6cOxc2boQnnoDHHtM20YGi5C8iIeezz8zU0IsvNp8KOnSwO6LIExEDviISWbp1MyuE\ne/Uyp4dt2GB3RNIU9fxFJGBWr4acHCgoMKUglYH8Q2UfEQl5n38OmZlwwQXw29+qDOQPEVH20Tx/\nkcjWqRP8/vfQp49ZGFZdbXdE4Uvz/EUkLC1bBqNHQ1ERjBmjdQG+UtlHRMLOp5/C8OHQs6fZLVRr\nArwXEWUfEXGWuDh47z3T6+/bF7ZvtzsiZ1LyF5Ggu+ACePllsy3EgAFw002wYAH87W92R+YcKvuI\niK1OnDCnhr36KqxbZ/YKGjkShgzRATJnopq/iESUv/0N3ngDXnzRJP4XX4TERLujCj1K/iISkerr\nYfFimDoVcnPh5z+Htm3tjip0aMBXRCJSixbw4IPw0Udmr6CePWH9erujiiwB7fnv2LGDOXPmcPz4\ncYYMGcKwYcMa31w9fxHxwNtvm+2ihw6F55/X+oCwKfscP36c+++/n6VLlza+uZK/iHjo0CEzK2jC\nBLNfkJMFreyTk5NDbGwsid8beamqqiI+Pp64uDiKi4ub/N0VK1YwYMAA7rrrrmYFKiLOdtFF8Lvf\nwaRJsGuX3dGEP496/hs2bCA6OpqRI0dSW1vbcL1Pnz7MmTOHLl26kJ6ezsaNG1mzZg01NTVMnDiR\nyy+/vOFnf/KTn7BixYrGN1fPX0S8VFxsNojbuBFat7Y7Gnv4I3e28uSHUlJScLvdja7V1dUBkJqa\nCkBaWhrV1dVkZWWRlZUFQGVlJaWlpViWxYgRI5oVqIgImIVhq1fD9OkwY4bd0YQvj5J/U7Zs2UL3\n7t0bHickJLB582aGDBnScK1///7079//rM9z6g51LpcLl8vla0gi4gBRUbBkidklND3dHBwf6Soq\nKvy++7HPyd9f/LU9qYg4x2WXwcKFkJUF27ZBTIzdEQXW9zvG06ZNa/Zz+jzPPykpiR07djQ83r59\nO3379vX6ebSfv4j4IiMDBg82U0Cdwpb9/N1uNxkZGU0O+F555ZUMHjyYjRs30r59e89vrgFfEWmG\nr7825wTfeCOMGmWmgjphDUDQpnpmZmbSr18/du7cSefOnVmyZAkAs2fPJi8vj0GDBjF27FivEv9J\n6vmLiK8uuAAqK6F7d3NIzDXXQGFh5E4F1UleIiLfY1lQUwOvvAKvvQbXXw8rVkB0tN2R+V/YrPA9\n482V/EUkAE6cMBvCtW4NixbZHY3/RcTGbir7iIi/tW4Nc+dCRQW8+abd0fiPyj4iIh7YvNlsBvfh\nh9Cpk93R+E9E9PxFRAKlb1+zIviBB8wZAfJPtid/lX1EJJAmT4YjR+C55+yOpPlU9hER8cKePfCj\nH0FZGfTubXc0zaeyj4iIB666yvT877nHLAyTEEj+KvuISDDce6/p9d9xB/z5z3ZH4xuVfUREfHDs\nGMycaY6C/MUvzL5ALVvaHZX3tMhLRMQHn3wCeXnmzeDFF80B8eFENX8RER/Ex5sFYLm5MHCgmRHk\ntKmgtid/1fxFxA4tWsCDD0JtLSxbZt4MQp1q/iIifvT00/CXv5jzgcOBav4iIn7w8cfmSMj//d/w\nOA9ANX8RET+IjzdnA9TU2B1J8Cj5i4jjRUXB7bfDW2/ZHUnwKPmLiGAWfy1bZncUwWN78tdsHxEJ\nBT/6ERw4AJ9+anckZ6bZPiIiATBmDFx9NUycaHckZ6cBXxERP3JS3V89fxGR7xw/DrGxZupnx452\nR3Nm6vmLiPhRmzbwb/8Gb79tdySBp+QvInIKp5R+VPYRETnFoUNwxRXw+edw8cV2R9M0lX1ERPzs\noosgJQVWr7Y7ksCyPflrnr+IhJpQXfClef4iIgH0f/8H3bvDl1/CeefZHc3pVPYREQmA2Fjo0QPW\nr7c7ksBR8hcRacLtt8PcubB/v92RBIaSv4hIE0aNMgu94uLMdg9ffml3RP6l5C8i0oQf/AAWLYJt\n2+DoUUhIgHHjzIEvkUDJX0TkLK680hzv+MkncOGF0KeP+VQQyrt/ekLJX0TEA7GxMHOmSfqdO0O/\nfpCZaQ6AD0dK/iIiXmjXDoqKYPdu8ykgLc0MDm/ZYndk3gl48j98+DBJSUmsWrUq0LcSEQmaiy6C\nxx83bwIDB8Kdd5pD4Kuq7I7MMwFP/rNmzeLuu+8O9G1ERGzRti3k58OuXTBiBOTkmO0h1q2DUF7D\n6lHyz8nJITY2lsTExEbXq6qqiI+PJy4ujuLi4tN+r6ysjISEBDp06OCfaEVEQlSbNpCbCzt2mBPB\nHnsMkpLMDqH19XZHdzqPtnfYsGED0dHRjBw5ktpTRjf69OnDnDlz6NKlC+np6WzcuJE1a9ZQU1PD\nxIkTeeGFFzh8+DAff/wxbdu25a233iIqKuqfN9f2DiISoerrYflyeOopM1V0yhS46y5o1ar5z+2P\n3Onx3j5ut5uMjIyG5F9XV4fL5WLr1q0AFBQUkJ6ezpAhQ0773ZdffpkOHTpw6623ntaAwsLChscu\nlwuXy+VrW0REQo5lwTvvmDeBL76AJ56AkSPNJwVPVVRUNNoAc9q0afYl//LychYvXszSpUsBmD9/\nPnv37mXGjBme31w9fxFxkKoq8ybwySdm1XBurhkz8FZEbOymLZ1FxClSU81A8JtvwrvvwtVXw6xZ\n5gAZT9iypfO5yj75+fkMHjy4ybLPGW+unr+IOFhtLfzyl1BWBg8/DAUFZh3Budja84+JiQHMjB+3\n201ZWRnJycnNCkZExEkSE+F3v4M//MEcG/n008G7t0fjzpmZmVRWVrJ//346d+7M9OnTyc7OZvbs\n2eTl5XHixAkKCgpo37691wEUFRVpoFdEHC0uzmwid67O/PcHfptDJ3mJiIQZDfiKiDiIzvAVEXGw\niOj5i4hI8Nme/FX2ERHxjMo+IiIOprKPiIj4RMlfRMSBbE/+qvmLiHhGNX8REQdTzV9ERHyi5C8i\n4kC2J3/V/EVEPKOav4iIg6nmLyIiPlHyFxFxICV/EREHUvIXEXEg25O/ZvuIiHhGs31ERBxMs31E\nRMQnSv4iIg6k5C8i4kBK/iIiDqTkLyLiQEr+IiIOZHvy1zx/ERHPaJ6/iIiDaZ6/iIj4RMlfRMSB\nlPxFRBxIyV9ExIGU/EVEHEjJX0TEgZT8RUQcKKDJv6KigpSUFMaMGUNlZWUgbxWSIn3xmtoX3iK5\nfZHcNn8JaPJv0aIF0dHRnHfeeVx99dWBvFVIivQ/QLUvvEVy+yK5bf7iUfLPyckhNjaWxMTERter\nqqqIj48nLi6O4uLi034vJSWFNWvWMGHCBJ555hn/RPw9nr7IZ/u5pv7Nk2unPj7T981lV/s8bWtz\nBat9drx2nj6ft21r6nok/W02dT2S2hcqucWj5J+dnc3atWtPuz5+/HgWLFhAeXk58+bNY9++fbz6\n6qs88sgjfPHFF0RFRQHQrl07Dh8+7LegTxXpL5CS/7l/Tsk/sv42m7oeSe0LldyC5aE9e/ZYPXr0\naHh88OBBq3fv3g2P8/PzrZUrVzb6ndLSUisvL88aNWqU9cEHH5z2nIC+9KUvfenLh6/maoWPtmzZ\nQvfu3RseJyQksHnzZoYMGdJw7Y477uCOO+4443NY2tRNRMQWmuopIuJAPif/pKQkduzY0fB4+/bt\n9O3b1y9BiYhIYPmc/GNiYgAz48ftdlNWVkZycrLfAhMRkcDxKPlnZmbSr18/du7cSefOnVmyZAkA\ns2fPJi8vj0GDBjF27Fjat28f0GBFRMQ/bD3JS0RE7BGSA7579uwhNzeXESNG2B2KXx07doxHH32U\nMWPGNLluItxF6ut20vLlyxk9ejQ5OTm8//77dofjdzt27GDMmDGMGjWK0tJSu8MJiMOHD5OUlMSq\nVavsDsXvKrzcTickk/9VV13FokWL7A7D7zZt2kRSUhIlJSUR+R9XpL5uJw0dOpSFCxfyq1/9qqH0\nGUm6d+9OSUkJJSUlvPHGG3aHExCzZs3i7rvvtjuMgPB2O52AJn9ft4UIJ960sba2lm7dugFw5MiR\noMfqi0h/DX1p38yZM8nLywtmmD7ztn0rVqxgwIAB3HXXXcEO1SfetK+srIyEhAQ6dOhgR6g+8aZ9\nXm+n0+xlYmdRVVVl1dTUNFoZbFmW1bt3b6uystJyu93WddddZ/31r3+1XnnlFWvChAnW3r17G35u\n+PDhgQzPL7xp4/r1662lS5dalmVZo0ePtiNcr3nTvpPC4XU7ydP27du3z6qvr7cmTpxolZeX2xSt\n93x5/SzLsjIyMoIZps+8ad/UqVOtCRMmWGlpadbQoUOt+vp6m6L2nC+vX11dnTVq1KhzPrfPK3w9\nkZKSgtvtbnStrq4OgNTUVADS0tKorq4mKyuLrKwsAA4cOMCUKVPYtm0bM2fOZNKkSYEMs1m8aeOg\nQYOYOnUqmzZtYtiwYcEO1SfetO/GG28Mm9ftJE/bt3nzZnbv3s369es5dOgQu3btCovevzevX3R0\nNKWlpViWFTbjNt6078knnwTg5ZdfpkOHDg17j4Uyb9p3/Phx1q1bxzfffMOYMWPO+dwBTf5N8WRb\niHbt2jF//vxgh+Y3Z2tjoHY3DaaztS+cX7eTztS+GTNmkJ+fb2Nk/nG29vXv39/GyPzjXDnm/vvv\ntys0vzjb63e27XS+LyQHfEVEJLCCnvydsC1EpLdR7Qtval9481f7gp78nbAtRKS3Ue0Lb2pfePNb\n+/w9On2qn/70p1bHjh2tNm3aWJ06dbJeeukly7Isq6KiwurevbvVrVs3a86cOYEMIeAivY1qn9oX\nytQ+39un7R1ERBxIA74iIg6k5C8i4kBK/iIiDqTkLyLiQEr+IiIOpOQvIuJASv4iIg6k5C8i4kBK\n/iIiDvT/1Kwh0Bt+XUIAAAAASUVORK5CYII=\n"
      }
     ],
     "prompt_number": 43
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "d = genfromtxt('fires.txt')\n",
      "x = powerlaw.Fit(d, discrete=False)\n",
      "x.sigma"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Calculating best minimal value for power law fit\n"
       ]
      },
      {
       "output_type": "pyout",
       "prompt_number": 29,
       "text": [
        "0.050979497997097889"
       ]
      }
     ],
     "prompt_number": 29
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print fit.power_law.alpha\n",
      "print fit.distribution_compare('power_law', 'lognormal', normalized_ratio=True)\n",
      "print fit.distribution_compare('power_law', 'exponential', normalized_ratio=True)\n",
      "print fit.loglikelihood_ratio('power_law', 'stretched_exponential', normalized_ratio=True)\n",
      "#print fit.loglikelihood_ratio('power_law', 'truncated_power_law')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1.63979125376\n",
        "(-7.1333027650672678, 9.7988686627588056e-13)"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "(11.477220800562032, 1.7171157386991705e-30)"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "(-7.1103191942301169, 1.1577488148942332e-12)"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ax = fit.plot_ccdf()\n",
      "fit.power_law.plot_ccdf(fit.data, ax, color='r')\n",
      "fit.lognormal.plot_ccdf(fit.data, ax, color='g')\n",
      "print fit.distribution_compare('power_law', 'lognormal')\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(0.92801788116372386, 0.42584569440313169)\n"
       ]
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEHCAYAAABGNUbLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlYVdX+x/H3AcQ5HHPWa+pNSUsrxAlERbEccMKhQgVH\nUH8OZaVWaDZpXpPIcNacFYdETVQyQFFRI8u8DteBciqnJHNWzu+Pc69EioEezoZzPq/n4amz2ez9\naT30dbn22muZzGazGRERcShORgcQERHbU/EXEXFAKv4iIg5IxV9ExAGp+IuIOCAVfxERB6TiLyLi\ngFT8RUQckEtOXfjGjRuMGjWKa9eu4e/vT+vWrXPqViIikk051vNPTEzEw8ODyMhIVq1alVO3ERGR\nh5Ct4h8cHEyZMmWoU6dOhuMJCQnUqlWLGjVqEBERAcC+ffuoVq0aANeuXbNSXBERsYZsFf+goCBi\nYmLuOT506FCmT59ObGwsU6dO5fz58zz99NMcO3YMgEKFClknrYiIWEW2xvy9vLxISUnJcCw1NRUA\nb29vAFq1akVSUhK+vr6MGTOGxMREOnXqdN/rmUymh4gsIiKPuibnI4/57969m5o1a9797O7uzs6d\nO8mfPz+TJk0iIiICPz+/TH/ebDZb9SssLMyq52Z2TlaPDx44mIa9i/Fsv1IcOnw0w/cz+3d7bYsH\nfbbntsjKMbWF7dsiu9fLTW1hDYZP9Rw7dixxcXFWu56Pj49Vz83snKwe79ytM1vCT/HMpQo0nVYT\n8rlkeq615ba2eNBne26LrBxTW9z/c062RXavnRvaIi4ujrFjx/5tjiwxZ9Px48fNtWvXvvv50qVL\n5rp16979PHjwYPO6deuydK2HuH3elZZmnhwUZC450sk8/KPXzGlpaRm+HRYWZkyuXEhtkU5tkU5t\nkc4atfORe/5ubm6AZcZPSkoKmzdvxtPT81Eva39MJobPmcOi4p+zOeUTWr7hzZWbV+5+O6d7e3mJ\n2iKd2iKd2sLKsvMnRffu3c3lypUzu7q6mitWrGieM2eO2Ww2m+Pi4sw1a9Y0V6tWzRweHp7l6wHm\nsLAw8zfffJOtP7HyukNrdpi7dCxs/serj5sPnj1kdBwRySO++eYbc1hYmFV6/iaz2UpPDx6CyWSy\n2sOLvObC/jOED2hMRJPTTO88n64eXY2OJCJ5hDVqp4q/ga5fus58X3/CmsbR9bn+/Kv7J7g45diK\nGyJiJ6xRO+1utk9eUqBYAfrtimHqv18leeNMvP7ViLNXzhodS0RyKWvO9lHPP5dIGL6MVfuCWNK4\nMNH91uFZUQ/NReT+NOxjZw4u3MPOD/wY1vEG77ebSKhniN6CFpF7qPjboYv7TrGnTSsGdjiNZ8MX\nmd1lJoXyaW0kEUmnMX87VKJOBZr/uJv5G5tzefVGPCI8OHrxqNGxRCQX0Ji/AzCnmUl48X0SL03i\nX22dmBfwBe2ebGd0LBHJBTTs4wDiR6zh6oreBAc70cc7hHE+43B2cjY6logYSMXfQcR//iOF3mzD\n/w0wU7TOkyzusoRShUoZHUtEDKLi70CSvrrAtY6dWdbpBBvq32Zlt1U8V/45o2OJiAH0wNeBeL5Y\nkop7N9Mo5kXeWnKb1l+05Iu9XxgdS0RsSA98HdjFi7DQewbPXxxFz8EFaV23A5P9JuPq7Gp0NBGx\nEQ37OKi0NFjQfysNF3dhyGuluPKP4qzouoKyRcoaHU1EbMAuhn0k+5ycoNcsL5LH7mLSRy547b7C\n8zOeZ8eJHUZHE5E8Qj3/PG7xzCsUGdKba433MaTVeca3eJ/+z/XXshAidswuev564PtoXupXmMc2\nLCcl6WViIvPxafxE+q3tx/Xb142OJiJWpge+co+ff4aZL65mwPF+DHmnOqdKpLGy60oquVUyOpqI\nWJld9PzFOipXhnE/dGRZwDd8+s6vtDtSmPqz6hOfEm90NBHJhdTztzNmM0wYeZ5mU7vwS+vrDGh0\nnDe9RjHUc6ieA4jYCU31lEwtmneL6yHD8KgSQ6+Q/LhXqsfMdloeWsQeaNhHMvVy73w8t2Mqiy+M\nYuPb5zCdOkOj2Y04/ttxo6OJSC6g4m/H6taF1w/15d0a0UweeYiuxyvSYHYDNh3dZHQ0ETGY4cVf\nUz1zVokS8MnOhkS9tpuWky8wdUtNeq/uxYdbP9SQm0geo6me8lC+33WDH5uG8vwTifQMKUClx6sz\n138uRfMXNTqaiGSDxvwlW56pn5/KG2cx5+f/Y1PYaYpfuEqD2Q04fOGw0dFExMZU/B2Ml7eJlxND\n6euygsmv72XoJXeazGnC+sPrjY4mIjakYR8HdeYMDOt8gnH7OnKybQl6Pf8jg+oPZlSTUXofQCSX\n0zx/eSR37sC0T65R+q3+uD/xHX0G5aPyf58DFHEtYnQ8EcmExvzlkTg7w6DXClL72/msOduHdWPO\nUOTCFb0PIOIAVPwF96dMBO8bzvjyC5k09FuCLv2ThrMb8vWxr42OJiI5RMM+ctedOzC6x3H6ru/A\n0W7l6O2+lzebvKl1gURyGY35i9WlpcEHY67w9JRg/uF+iF5Bt3m64nNMazONgvkKGh1PRLCTMX+9\n4Zu7ODnBWx8WpvDapcQc7876ty5w/cwpvOd5c/L3k0bHE3FoesNXbCI1FSY228Br+3vy+XgvPndJ\nYnmX5TSu3NjoaCIOTcM+kuNu3YK3uv2H0E3+7A+uSlClPbzb7F0GPD/A6GgiDkvFX2zizh1o7/M7\n43/qScEaJ+jc/jLe1Zrz6Quf4ursanQ8EYdjF2P+kvs5O8OitY8xtcUqNu1tT/zHVzhz4gAt5rfg\n1z9+NTqeiDwEFX/JkmLFYPZcJ8pPC+PVPz4nasx+WqSWxGOmB3tO7zE6nohkk4q/ZEtAALSM8Kc5\nCYyY8G/CT9bmhUUvsOD7BUZHE5Fs0Ji/PJTkZBjQ7RILCOR29dN0aHWR9u4dmdhyIi5OLkbHE7Fr\nGvMXwzz7LKxPLEZwyTUk/diOHeE32f+f7bRe2JoLVy8YHU9E/oaKvzy0xx+HbdudODdoLKOuRhI9\n9ijPnnPBY6YH+37dZ3Q8EXkADfuIVcyaBYvePshqU0fWda3A8PLfM6v9LPxr+hsdTcTuaJ6/5Cob\nN8LHb//O2OO9MD1zlG6tzxNSfxCjvUZrYTgRK8r1Y/7Hjx+nb9++BAQE5ORtJJfw84OY7Y+xddhK\n4rZ3ZfPHaazatYgeK3tw9dZVo+OJyJ/YpOcfEBBAVFTUvTdXz99unToFX/ZfT7vNvXml3xNcrn2H\ntS9/ScXHKhodTSTPs1nPPzg4mDJlylCnTp0MxxMSEqhVqxY1atQgIiLikYKIfalQAQatb4Prlu0s\nXnwZ9ygnnov0ZMeJHUZHExGyWPyDgoKIiYm55/jQoUOZPn06sbGxTJ06lfPnz7NgwQKGDx/O6dOn\nrR5W8p6yTWpQ8eck3vu9MqPmFqHFjHZEbv/C6FgiDi9Lxd/Ly4vixYtnOJaamgqAt7c3VapUoVWr\nViQlJREYGMgnn3xC+fLluXjxIgMHDmTv3r1MmDDB+uklbyhalKq7owjpEcSWL0y8vnIMA1e9xp20\nO0YnE3FYD/0q5u7du6lZs+bdz+7u7uzcuZM2bdrcPVaiRAmmTZv2wOv8eWMCHx8ffHx8HjaS5GYm\nE/nD3sSzfl32BwbS/Pd1+Jzdz9reSyhWoJjR6URytbi4OKtvemX4e/jW2pVG8gbTC62ptHMnW306\nEPRLCk+e9WTroLX8s+Q/jY4mkmv9tWM8bty4R77mQ0/19PDw4ODBg3c/79+/nwYNGmT7OtrG0fGY\nqlej3KGdrOQZQtZfx+Ozxmw6usnoWCK5niHbOKakpNCuXTv27Ut/bb9evXqEh4dTuXJlWrduzbZt\n2yhVqlTWb66pno7NbObXN/7FzqgP6fUKjG3zNkM9h+qFMJG/YbOpnj169KBRo0YcPnyYSpUqMXfu\nXACmTJnCgAED8PX1JTQ0NFuFXwSTiTITX8N74jK2fO5E+MqP6bWyLzdu3zA6mYjdM3x5h7CwMD3o\nFU4lpnCpU3v6NL/Edc8KbOz7JWWKlDE6lkiu8r8Hv+PGjdPaPmJHrl7lYLN+hBfewnIvZxZ0WMuL\n9eoZnUok18n1a/uIZEuhQtTcuZAwzzf4ZM0fdFzajOJNoligTcJErM7w4q/ZPpKByUTZD4fR85PV\n7FzvQuEW/Rm4NIyFi9JISzM6nIixDJntkxM07CMP9PPP/Nq9Hf4NTnP0ShPK7l5I9MrCVK1qdDAR\nY2nYR+xb5cqU+Xon8Rf8ePFWHL+3a8CzPidZuhTUZxB5NOr5S+5nNmP+7DM+XjeKyY0LYor+iqbV\nPZgxAx57zOhwIrZnFz1/jfnL3zKZMA0ZwutjvmJazB1ut2vGryWX0bw5TJoEy5ej5wHiEDTmL47r\n5En29m6Nv+dxni8+jCpn3iNquYmgIBg3DvRysDgCu+j5i2RLxYrUXbeHpLPtOfXzJ5yu3ZaEHdeI\nioJ33oEbejlYJEtU/CXvKVCAsjMWE1drAs6bv6bbkmeYMusMW7dCz57w0096ICzydwwv/hrzl4di\nMlEgZAgLB32Nf/wv9IutxYdz95AvHzRoAM88A99/b3RIEevSmL/In50+zcpBzRhYO4VpHWbR6dlA\nZs2CsWNh/nxo0cLogCLWZY3aqeIv9uHGDZJffQn/wtEM9BzM6I6TWbfORGgo1KkDI0aAr6/RIUWs\nQ8Vf5C9OT5+E//ejeLJmY2aFxmC6U4BFi2DMGJg3D/z8jE4o8uhU/EXu4+r2eIKmv8DPVUvw5fBd\nlHErT3w8dOsGkydDly7g6mp0SpGHp6meIvdRqFFTlk44hl+KC54TqvP94a00bQoLFkBkpOVh8J92\nIBVxSIYXf832kZxgKluWsTP/w0e3ffCd60P0pghatoT4eBg4EJo3hz/+MDqlSPZoto9INuya/S4d\nD45jaNXujAxZiMlkYuBAy1TQLVugYEGjE4pkj8b8RbLo5I6NtF/Snqcfq8H0t5LI51qYDh0gIQH2\n7IHq1Y1OKJJ1GvMXyaKKDf3YOuYIl38/R4vRFbhw8gDR0ZYlIZ5+GkJDjU4oYlsq/uIwCpepRNTk\nkzQt/BSe4XX49zfLGTECfvkFNmyA11/XcwBxHCr+4lCcXPLx/vhExrqH4hPTndjIkTz2GMTGwuHD\n0LQppKYanVIk52nMXxxWQvx8um4I5t1bTej/4SbM+Vzp0wd27LA8CC5XzuiEIvdnF2P+muopRvFu\n2pOtA3fxL9c9vDqgKmmnTjB7NnTsCMOGGZ1O5F6a6iliRRevnKfzpPo8lnKGRb2juVqrJU89Zdkl\nLDAQnAzvIolkpKmeIlZy885NBk5ry3cHvmHtk2M5XX80ffqaqFsX5s4FFxejE4qks4thH5HcwNXZ\nldmhG+nWchgNT47DZY4/ibHXOH8eWraEr76Cs2eNTiliPer5i/zFiuSFhKzuy6zkirSN+JppG6qw\nciXs3g3R0dCsmdEJxdFp2Eckh+w6mUTHOa0YkZjGiNe/xNSiBaNHw4EDMG0alCljdEJxZBr2Eckh\n9St6smPoPr7wLc3AGe249fEERgw3U6GCZSkIbREpeZ16/iIPcPnGZbov9Ofmd3uIuuhLsZkLmPBZ\nYRYtgmXLoFYtoxOKI1LPXySHFc1flDVBm3BvHUijqls41uJZhrU7Qvv2lreBV6wwOqHIw1HPXySL\nPkuK4P1Nb7FyhTONPlzI7tIv0qWLZYewiRONTieOxC56/nrDV/KKwZ5DmNN9KR26mVky4RU8Ysbz\n/XdpLFoEPXpYHgaL5CS94StioH2/7qPdwhcJ2mvind/rcm7SAmYsc+PDD+Gzz6B1a60LJDlLUz1F\nDPLLH7/gv7g9NY6nMuvLNPKvWMPCZHfWr4dt22D/fnBzMzql2Cu7GPYRyYvKFilLXFA8N599Bt9e\nJs639iKw0EqWLoU2bcDX17I20JUrRicVuT/1/EUeQZo5jbe2vEVU8kK+mn+bGm17cfOd94jZ7Myk\nSZZzZsyAmjWNzSn2RcM+IrnErORZvBU7mqgdlfC6WgoWL+Z64ZKMHw/ffgsxMUYnFHuiYR+RXKLv\ns32Z33khnRudYPGz+cDDgwIH9zJmjGUW0LZtRicUyUg9fxEr2vfrPtouaUv/fA0ZPTYW05RwIn9/\nmchISE7W0tBiHer5i+QydcrUYUefHax0OUzfj725NfYdQg4Nw9V0i82bjU4nkk49f5Ec8MfNP+i+\nojs3rl9hxZf5uHn0Bu2vLWfHMS0HKo9OPX+RXKqIaxG+7P4lT5Z9isatT3Ol0zNEpTxPUkSS0dFE\nABV/kRzj4uRCxAsR9KnXh8bFV7Jt7DCqDWvHhBozOXrU6HTi6HJ02GfNmjWsX7+e27dvM3DgQOrX\nr5/x5hr2EQex+sBq+q/rz8y642nQ/VP+U6YJTZIjMBXIb3Q0yYPyzDz/s2fPEhYWRmRkZMabq/iL\nA9l1ahcdlnagT5XhNO69k9I3TjG3zQoGvleR2rWNTid5ic3G/IODgylTpgx16tTJcDwhIYFatWpR\no0YNIiIiMv35CRMmMGDAgEcKKpLX1a9Qn8TgRFb8OoeNKypTZUR7PtxSnw9eSODOHaPTiaPJUs9/\n69atFClShJ49e7Jv3767x+vVq0d4eDhVqlTBz8+Pbdu2sWHDBpKTkxk5ciTlypXjjTfewM/PjxYt\nWtx7c/X8xQH9du03Oi3vhFt+Nxa5BXO9fT9O9BxD3VlDwGQyOp7kAdaonVl65cTLy4uUlJQMx1JT\nUwHw9vYGoFWrViQlJREYGEhgYCAAn376KVu2bOHy5cscOXLkvr3/P69N7ePjg4+Pz0P8Z4jkHcUL\nFmfjKxvpG90Xn5PjGRmxlqeG9Gftyt0cGjGdQSMLUbCg0SklN4mLi7P6vidZHvNPSUmhXbt2d3v+\nsbGxzJ49myVLlgAwbdo0Tp06xfjx47N+c/X8xYGZzWbejX+Xed/PY22HFTze/xOu7NpPsNsq5m+t\nSqVKRieU3MpmPX8RsT6TyUSYTxhVi1el+fIXWBq5hOZrf2TtmAYMa7yAUi+1okkTaNQISpQwOq3Y\nm4ee5+/h4cHBgwfvft6/fz8NGjTI9nW0jaM4up7P9GR5wHJ6rHqJeU3dKBi9nM+u9MZ3z0eETzHz\nj3/AlClGp5TcwJBtHP867APpD3wrV65M69at2bZtG6VKlcr6zTXsI3LXgXMHaLO4DS8//TLvVu+P\nqUsXqFSJ5CFz6dG/KIcOGZ1QcgubTfXs0aMHjRo14vDhw1SqVIm5c+cCMGXKFAYMGICvry+hoaHZ\nKvwiklGt0rXY0WcHm45uInD3KG5+EwvFi1MvtAEuxw5z86bRCcWeGL6wW1hYmGb5iPzJ1VtXeXnV\ny6ReT2VVt1UUWxDF+YFjuPDRLJ4c2d7oeGKg/836GTduXN54wzfTm2vYR+S+7qTdYdjGYXxz/Bs2\nvLyB/VNO8fS7AZz3D6ZOVBgmZy3L5cjyzPIOmd5cxV8kU2azmck7JjMlaQrreqzj5p6yOHUPoGDp\norh/twiKFTM6ohjELpZ01mwfkfszmUy82uhVJrWcRMsFLblU6wf+ceRrEn+tztl/ePD9oh9R38mx\nGDLbJyeo5y+SNQk/JRAQFcBE34k0L9mL5BEL8Vo9nHGlp9JoSle6dTM6odiShn1EHMiBcwd4cfGL\nBNcN5i3vt+C7vdxo24kF17qwt+uHfDDRBTc3o1OKLaj4iziYM5fP0HZJW+qVrUdkm0jyXfqdW11f\n4j8H79DpxlL6jSrFoEFQoIDRSSUnacxfxMGUK1qO+N7xnLp8ivZL23O5iCv5Nn2Fe08Pfsj/PCe+\n/JaGDeHECaOTSk7QmL+Ig7t15xahX4Xy7elvWf/SesoVLQcrV2IOCSGm+cf03dqL5GQoo/3i7ZJd\n9PxFJPvyOedjRtsZdKrViYazG3Lg3AHo3BlTXBwvfPcBcwoPZt0qvRIsmVPPXySPm//9fEZuHklU\nQBTeVbwhNZVf/Hpycu8FqiVHUdy9nNERxcrsouevMX+RR9PzmZ4s7LiQLsu7sOzHZeDmRtntqzn3\nrB9pz3vA9u1GRxQr0Zi/iNzj+1++p+2Stgz1HMqrDV/l1CkTo575igUuQTB2LAwcqG0i7YSmeopI\nBidST/Di4hdpWqUpHzcPp2QJZy4kHaHgSx3BwwM+/1zzQO2AXQz7iIj1VHKrxNagrfz73L/p8WVn\nnvW8StzJ6rBjB1y5Al5e8PPPRseUXEDFX8TOFCtQjJhXYiiavyg/N2vO6k3noEgRWLoUunUDT0/4\n5hujY4rBVPxF7JCrsyvzO8ynZbUWzHNuxJGLRyzj/a+9BgsXwksvweTJaGU4x2V48ddsH5GcYTKZ\nmNnjfQrtfY2GM73YdWqX5RstWsDOnbB4seUPgStXjA0qWabZPiKSZW+8AUec15JQPJh5/vNo8882\nlm9cuwahofDtt7B6NVSrZmxQyTI98BWRv9WnDyTMbMe4J9fRJ7oPs5JnWb5RsCDMmWOZAtqoEWzY\nYGxQsSn1/EUcwJYtlme9k+YeZtzR1vSu25u3vd/G9L95/4mJ0LUrhITA6NHgpH5hbqZ5/iKSZfHx\n0KULBAT9wtbKbWhY5Tk+b/M5Lk4ulhNOn4aAAChdGubPh8ceMzawZErDPiKSZU2bwu7dUPBOWU69\nF0d0/E80jezIHzf++8C3fHnLFNDy5aF+fThwwNjAkqPU8xdxQJcvw5x5t3h7d1/SShwi3HMtfXqU\nTj9h7lzLk+Lp06FjR+OCyn3ZRc9fUz1FbK9oURg6JB+/zZ1Hm1otGLCrMeu3H0s/ISgIvvoKhg2D\nMWPgzh3jwspdmuopIlbV4f3P2XT9fbaFruXZcs+mf+PcOcuTYldXy3sBJUoYF1Lusouev4gYb96g\nUJxjPqPlF63ZdHRT+jdKl4ZNm6B2bXj+efj+e+NCilWp+IsIxYrBp6EdubNoNV0WB/LF3vnp33Rx\ngUmT4IMPwNcXFi0yLqhYjYZ9ROSu/fuhy8AD/OT1AsO9BvBe6zfT3wUA2LfP8gC4bVv4+GPIl8+4\nsA5M8/xFxOpu3oSR757m89QXaFGjCWtDPyWfi3P6Cb/9Bq+8YlkTaNky7RJvAI35i4jVubpC+Hvl\n+SoggcRDBynevxvjP7zOL7/894TixWHtWvD2tjwHSEoyNK88HPX8RSRT12/doO2cnuxP+ZVrc9bQ\nvLEb/fpBq1bg7AxER0PfvvD++9Cvn9FxHYaGfUQkx6WZ0xgaM5T441sJNG1g+ZxynD0LwcGWr0pX\nD1meAzRpAhERkD+/0ZHtnl0M++glL5HczcnkxKetP6V7na5E3mzMopjDrFkDZ89C3brQZsST7J+T\nBBcvWoaCTp40OrLd0kteImKI2cmzeeubt4juHo1HBQ+uXoVp0yAyEvZ+Z6bwZxPg009hyRLLYkKS\nIzTsIyI2F30omj7RfVjYcSF+1f0ACAwENzf47DMsL4X17AmjRsH//Z9l+0ixKhV/ETFE4s+JdFre\nicmtJvPy0y9z6RLUqWPZG6ZlS+D4cejUCdzdYeZMKFTI6Mh2xS7G/EUk72lcuTFbem5h1NejmLxj\nMsWKWQp/cLDlNQCqVrVsEOPkZNkl7Nixv72m2JaKv4g8lKcef4rE4ERmJc/i9c2v08I3jRYtLM8A\nAEtvf/58yz6SDRvCxo2G5pWMNOwjIo/k4rWLtF3cluolquNxejY//pCP6dP/clJCAnTvDoMHW54F\n6DnAI9Gwj4gYrkTBEsT2jOW3678x74Y/P525cu9J3t6WbcTWroXOneH3320fVDJQ8ReRR1YoXyFW\nd1tNBbeybK3anPNXz997UoUKEBcHjz8Onp5w8KDNc0o6FX8RsQoXJxemvTAbjregyZwmpFxKufek\n/PktDwVee83yt4Evv7R5TrFQ8RcRqylTxsTtjR/QtWooTeY04Ydff7j/iX36wLp1lvcAtE2kIVT8\nRcRqnJ1hyhSYFvR/dC/xL3zn+xKfEn//k+vXhz17YPt2y/4AFy/aNqyDU/EXEasKCbE8140K64bP\nhcUERAWw6sCq+5/8+OOweTPUqgUeHtom0oZydKrnwYMHCQ8P5+bNm7Rp04ZOnTplvLmmeorYrQsX\nLMs+nCGZM83aMrbZOwx8fmDmP7B4MQwdCuHh8NJLtguaB+WZ5R1u3rxJr169WLJkScabq/iL2LW0\nNPjwQ5gy/yiuwX709wzknabvZNwa8s9++MGyPHT79jBxoraJzITN5vkHBwdTpkwZ6tSpk+F4QkIC\ntWrVokaNGkRERNz3Z6Ojo2nWrBldu3Z9pKAikvc4OVme5y6LrMbtGYlMj1/DwHUh3EnL5AHv009b\n3gc4eNCySNDZs7YN7ECy1PPfunUrRYoUoWfPnuzbt+/u8Xr16hEeHk6VKlXw8/Nj27ZtbNiwgeTk\nZEaOHEn58uXvntu+fXuio6Mz3lw9fxGHcfo0dHn5dw7V7USjZ92I6rGIAi4F7n/ynTsQFmZZHmLF\nCsvDYbnLGrXTJSsneXl5kZKSkuFYamoqAN7e3gC0atWKpKQkAgMDCQwMBCA+Pp5Vq1ZhNpsJCAi4\n77X/vDGBj48PPj4+2fxPEJG8oHx5SNj8GG++tZ7PY3vS/LofGwKjcSvgdu/Jzs7w3nuWPYLbtoUP\nPrBsF+mg4uLirL7pVZbH/FNSUmjXrt3dnn9sbCyzZ8++O44/bdo0Tp06xfjx47N+c/X8RRxSo8Zp\nlA8eytHb24h5OYYyRcpkfvLBg5bnAN7elo1itE2k1vYRkbypVEknXin5KR2e7ECTuU04/tvxzE+u\nWRN27YLz5y27g2mbSKt46OLv4eHBwT+tzbF//34aNGiQ7etoD18Rx1OiBPz2m4kwnzCGeQ7Da64X\nP579MfMfKFrUMvbv728Z/09IsF3YXMSae/g+dPF3c7OM0yUkJJCSksLmzZvx9PTM9nXGjh2rcX4R\nB1OyZPoLvYPqD+Ljlh/TYn4Ltp/YnvkPmUyW5aDnzoWAAMsQkIMNG/v4+Ni2+Pfo0YNGjRpx+PBh\nKlWqxNwhgBMSAAAJJUlEQVS5cwGYMmUKAwYMwNfXl9DQUEqVKmWVUCJi30qUsLwE9j896vTgiw5f\n4L/Unw3/2fDgH/bzg507LVuH9ewJV6/mbFg7ZfhmLmFhYZrlI+JgIiMtKznc3fXrv3ac2EGHZR2Y\n4jeFHnV6PPgiV69C//6wfz+sWmXZOtLO/W/Wz7hx4/LGG76Z3lyzfUQc0vLlEBVl+fqrH8/+SOuF\nrXmzyZsMrj/4wRcymy3DPx98AAsWQKtWORM4l9FsHxHJk0qUyHwRz9qP12Zr0FbCk8IZGzf2wUXO\nZLKsB7R8OfTubVlLQh3KLDG8+Gu2j4jjKVky45j/X1UtXpVtQdtYc2gNQzYMIc2c9uALNm1qmQ66\nZg106QKXL1s3cC5hzdk+GvYREZs7fdqyivPXX1te4s1M6vVU2i9tT4WiFZjXYR6uzq4PvvCNGzBk\nCPz8M8TEWDd0LpJnVvXM9OYq/iIOa9kyGDwYRo+GYcMsIzj3c+3WNbqt6MattFusCFhBYdfCf3/x\n8+fBjmcfqviLSJ527Bh0727Z02XevMzr9e202/SN7svhC4dZ99I6ShQsYdOcuY1dPPDVmL+I43ri\nCdi2DdzdoW5dyKwUuDi5MMd/Dg0rNaTpvKacvnzapjlzC435i4jdiYmBoCDL1P233waX+6w5bDab\nmZA4gRnfzmBT4Caql6hu+6C5gIZ9RMSunDlj2frx5k3Lro4VK97/vJnfziQsLoz1L62nXrl6tg2Z\nC9jFsI+IyP+UKwcbN0Lr1vDcc/CX/Z/u6vdcPyJeiMBvoR8JPznmIm+PyvDirzF/EfkzZ2fLDKBV\nqyyzNocOtczg/KvO7p1Z0nkJXZZ3IfpQJn9K2BmN+YuIQ/jtN+jTB1JSYOlS+Oc/7z1n96ndtFvS\njgm+E+hVt5fNMxpBwz4iYteKF4eVK6FfP2jc2LJ8z195VPAgrncc78S9w+Qdk20fMo9Sz19E8oQf\nfoBu3Sx7uUydCkWKZPz+idQTtFrYio41O/J+8/cxZfbWmB1Qz19EHMbTT8OePZYpoM89B999l/H7\nldwqsTVoK7HHYhn01SBjQuYhhhd/PfAVkawqXBhmz4awMMvqzRERGRfxLFWoFF/3/Jr2T7Y3LmQO\n0gNfEXF4R45YloaoUMGyqVfJkkYnsh0N+4iIw6peHbZvt/yzXj3YutXoRHmLev4ikuetX2+ZEjpo\nkOUdAWdnoxPlLC3vICLyX6dOwSuvWP590SIoX97YPDlJwz4iIv9VoQLExkKLFpa/BciDqecvInbn\n1i3Il8/oFDnHLnr+muopItZmr4VfUz1FRByYXfT8RUTE9lT8RUQckIq/iIgDUvEXEXFAKv4iIg5I\nxV9ExAGp+IuIOCAVfxERB2R48dcbviIiWaM3fEVEHJje8BURkYei4i8i4oBU/EVEHJCKv4iIA1Lx\nFxFxQCr+IiIOSMVfRMQBqfiLiDggFX8REQek4i8i4oBU/EVEHFCOF/8rV67g4eHB+vXrc/pWeZoW\nt0untkintkintrCuHC/+EydOpFu3bjl9mzxPv9jp1Bbp1Bbp1BbWlaXiHxwcTJkyZahTp06G4wkJ\nCdSqVYsaNWoQERFxz89t3rwZd3d3SpcubZ20WZCdX5CsnJvZOVk9/qDPOf3LrLbI/N6Pem522iIr\nx9QW9/+ck22R3WvbW1tkqfgHBQURExNzz/GhQ4cyffp0YmNjmTp1KufPn2fBggUMHz6c06dPEx8f\nz86dO1m8eDEzZ860yfLNKniZ3/tRz1Vb/P05ufF/8qzkeZRz82pbOHrxx5xFx48fN9euXfvu50uX\nLpnr1q179/OQIUPM69atu+/Pzps3z7x+/fp7jgP60pe+9KWvh/h6VC48pN27d1OzZs27n93d3dm5\ncydt2rS559xevXrd9xpmbeQiImIITfUUEXFAD138PTw8OHjw4N3P+/fvp0GDBlYJJSIiOeuhi7+b\nmxtgmfGTkpLC5s2b8fT0tFowERHJOVkq/j169KBRo0YcPnyYSpUqMXfuXACmTJnCgAED8PX1JTQ0\nlFKlSuVoWBERsQ6TWU9dRUQcTq574Hvjxg1GjBhBSEjIfd8tcCTHjx+nb9++BAQEGB3FcGvWrKF/\n//4EBweza9cuo+MY6uDBg4SEhNCnTx9WrVpldBxDafmYdHFxcXh5eRESEkJ8fPzfnp/rin9iYiIe\nHh5ERkY6/C921apVmTVrltExcgV/f39mzJjBRx99dHfY0VHVrFmTyMhIIiMjiYqKMjqOobR8TDon\nJyeKFClC/vz5eeKJJ/7+fBtkytbyEPv27aNatWoAXLt2zRbxbOphl8qwRw/TFhMmTGDAgAG2jGkT\n2W2L6OhomjVrRteuXW0dNUdlpx2MWD7G1rLTHl5eXmzYsIFhw4YxadKkv7/4I78mlgUJCQnm5OTk\nDG8Im81mc926dc3x8fHmlJQU85NPPmk+d+6cecuWLeYlS5aYzWazuX///raIZ1PZaYv/6dKli61j\n2kRW2+L8+fPmtLQ088iRI82xsbEGpc1ZD/N7YTabze3atbNlzByXnXYYM2aMediwYeZWrVqZ/f39\nzWlpaQalzjkP83uRmppq7tOnz99e+6Hf8M0OLy8vUlJSMhxLTU0FwNvbG4BWrVqRlJSEr68vY8aM\nITExkU6dOtkink1lpy0aNmzI6NGj2bt3LxMmTOCNN96wddwcldW22LlzJ8eOHWPLli1cvnyZI0eO\n2F3vPzu/F0WKFGHVqlWYzWa7ex6UnXZ47733APjiiy8oXbo0JpPJplltITvtcfPmTTZu3Mjt27cJ\nCQn522vbpPjfz4OWh8jSX1nsyIPaYtq0aQYms73M2mL8+PEMGTLEwGS296C2aNq0qYHJbOvvlpLJ\nbPkYe/Wg34uOHTtm+Tq57oGviIjkPMOKv5aHSKe2SKe2SKe2sFA7ZGSt9jCs+Gt5iHRqi3Rqi3Rq\nCwu1Q0ZWaw9rP52+n+7du5vLlStndnV1NVesWNE8Z84cs9lsNsfFxZlr1qxprlatmjk8PNwWUQyn\ntkintkintrBQO2SUk+2h5R1ERByQHviKiDggFX8REQek4i8i4oBU/EVEHJCKv4iIA1LxFxFxQCr+\nIiIOSMVfRMQBqfiLiDig/wcHVzCZswHMMAAAAABJRU5ErkJggg==\n"
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "a_param1 = fit.exponential.parameter1\n",
      "a_param2 = fit.exponential.parameter2\n",
      "a_like = sum(fit.exponential.loglikelihoods(fit.data))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%R -i data,xmin,a_param1 -o s_like\n",
      "source('pli-R-v0.0.3-2007-07-25/exp.R')\n",
      "#fit <- exp.fit(data, threshold=xmin)\n",
      "#s_param1 <- fit$rate\n",
      "s_like <- exp.loglike.tail(data,a_param1,xmin)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 13
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%R -i data,xmin,a_param1 -o s_param1,s_like\n",
      "source('pli-R-v0.0.3-2007-07-25/discexp.R')\n",
      "fit <- discexp.fit(data, threshold=xmin)\n",
      "data <- data[data>=xmin]\n",
      "s_param1 <- fit$lambda\n",
      "s_like <- discexp.loglike(data,a_param1,threshold=xmin)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 14
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print a_param1, s_param1\n",
      "print round(s_like[0],4), round(a_like,4)\n",
      "print round(s_like[0],4)==round(a_like,4)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "0.0183978968538 [ 0.01838506]\n",
        "-14778.8513 -14778.8513\n",
        "True\n"
       ]
      }
     ],
     "prompt_number": 15
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "a_param1 = fit.lognormal.parameter1\n",
      "a_param2 = fit.lognormal.parameter2\n",
      "a_like = sum(fit.lognormal.loglikelihoods(fit.data))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": "*"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%R -i data,xmin,a_param1,a_param2 -o s_param1,s_param2,s_like\n",
      "source('pli-R-v0.0.3-2007-07-25/lnorm.R')\n",
      "fit <- lnorm.fit(data, threshold=xmin)\n",
      "s_param1 <- fit$meanlog\n",
      "s_param2 <- fit$sdlog\n",
      "s_like <- lnorm.loglike.tail(data,a_param1,a_param2,xmin)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": "*"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%R -i data,xmin,a_param1,a_param2 -o s_param1,s_param2,s_like\n",
      "source('pli-R-v0.0.3-2007-07-25/disclnorm.R')\n",
      "fit <- fit.lnorm.disc(data, threshold=xmin)\n",
      "s_param1 <- fit$meanlog\n",
      "s_param2 <- fit$sdlog\n",
      "s_like <- lnorm.tail.disc.loglike(data,a_param1,a_param2,xmin)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": "*"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print s_param1, s_param2, a_param1, a_param2\n",
      "print round(s_like[0],4), round(a_like,4)\n",
      "print round(s_like[0],4)==round(a_like,4)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": "*"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "d = fit.data\n",
      "xmin = fit.xmin"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": "*"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%R -i d,a_param1,a_param2 -o s_cdf\n",
      "s_cdf = plnorm(d,a_param1, a_param2, lower.tail=FALSE, log.p=TRUE)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": "*"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "lognormal = powerlaw.Lognormal(xmin=xmin, parameters=(a_param1, a_param2))\n",
      "f = figure()\n",
      "ax = f.add_subplot(111)\n",
      "ax.scatter(d, s_cdf, color='r')\n",
      "ax.scatter(d, log((1-lognormal._cdf_base_function(d))), s=1)\n",
      "#ax.set_xlim(0, 2000)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": "*"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%R -i d,a_param1,a_param2 -o s_pdf\n",
      "s_pdf = dlnorm.tail.disc(d,a_param1, a_param2, threshold=xmin)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": "*"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "f = figure()\n",
      "ax = f.add_subplot(111)\n",
      "ax.set_yscale('log')\n",
      "ax.set_xscale('log')\n",
      "ax.scatter(d, s_pdf, color='r')\n",
      "ax.scatter(d, fit.lognormal.likelihoods(d), s=1)\n",
      "f = figure()\n",
      "ax = f.add_subplot(111)\n",
      "ax.plot(d, s_pdf-fit.lognormal.likelihoods(d))\n",
      "figure()\n",
      "hist(around(s_pdf-fit.lognormal.likelihoods(d), 5), normed=True)\n",
      "#ax.set_xlim(0, 2000)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": "*"
    }
   ],
   "metadata": {}
  }
 ]
}